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Abstract 


In  this  report  we  consider  an  electromagnetic  interrogation  technique  for  identifying  the 
dielectric  parameters  of  a  Debye  medium  in  two  and  three  dimensions.  These  parameters 
include  the  dielectric  permittivity,  the  conductivity  and  the  relaxation  time  of  the  Debye 
medium.  In  this  technique  a  travelling  acoustic  pressure  wave  that  is  generated  in  the  Debye 
medium  is  used  as  a  virtual  reflector  for  an  interrogating  microwave  electromagnetic  pulse 
that  is  generated  in  free  space  and  impinges  on  a  planar  interface  that  separates  air  and  the 
Debye  medium.  The  reflections  of  the  microwave  pulse  from  the  air-Debye  interface  and  from 
the  acoustic  pressure  wave  are  recorded  at  an  antenna  that  is  located  in  air.  These  reflections 
comprise  the  data  that  is  used  in  an  inverse  problem  to  estimate  the  dielectric  parameters 
of  the  Debye  medium.  We  assume  that  the  dielectric  parameters  of  the  Debye  medium 
are  locally  pressure  dependent.  A  model  for  acoustic  pressure  dependence  of  the  material 
constitutive  parameters  in  Maxwell’s  equations  is  presented.  As  a  first  approximation,  we 
assume  that  the  Debye  dielectric  parameters  are  affine  functions  of  pressure.  We  present  a 
time  domain  formulation  that  is  solved  using  finite  differences  in  time  and  in  space  using  the 
finite  difference  time  domain  method  (FDTD).  Perfectly  matched  layer  (PML)  absorbing 
boundary  conditions  are  used  to  absorb  outgoing  waves  at  the  finite  boundaries  of  the 
computational  domain  and  prevent  excessive  spurious  reflections  from  reentering  the  domain 
and  contaminating  the  data  that  is  collected  at  the  antenna  placed  in  air. 

Using  the  method  of  least  squares  for  the  parameter  identification  problem  we  solve  an 
inverse  problem  by  using  two  different  algorithms;  the  gradient  based  Levenberg-Marquardt 
method,  and  the  gradient  free,  simplex  based  Nelder-Mead  method.  We  solve  inverse  prob¬ 
lems  to  construct  estimates  for  two  or  more  dielectric  parameters.  Finally  we  use  statistical 
error  analysis  to  construct  confidence  intervals  for  all  the  presented  estimates.  These  con¬ 
fidence  intervals  are  a  probabilistic  statement  about  the  procedure  that  we  have  used  to 
construct  estimates  of  the  various  parameters.  Thus  we  naturally  combine  the  deterministic 
nature  of  our  problem  with  uncertainty  aspects  of  estimates. 
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1  Introduction 


The  study  of  transient  electromagetic  waves  in  lossy  dispersive  dielectrics  is  of  great 
interest  due  to  the  numerous  applications  of  this  subject  to  many  different  problems 
[AMP94,  APM89].  These  include  microwave  imaging  for  medical  applications  in 
which  one  seeks  to  investigate  the  internal  structure  of  an  object  (the  human  body) 
by  means  of  electromagnetic  fields  at  microwave  frequencies  (300  MHz  -  30  GHz).  This 
is  done  for  example,  to  detect  cancer  or  other  anomalies  by  studying  changes  in  the 
electrical  properties  of  tissues  [FMS03].  These  include  the  dielectric  properties  such 
as  permittivity,  the  conductivity  and  the  relaxation  time  for  the  Debye  medium.  One 
generates  the  electrical  property  distributions  in  the  body  ( Microwave  maps )  with  the 
hope  that  such  properties  of  different  bodily  tissues  are  related  to  their  physiological 
state.  With  such  noninvasive  interrogating  techniques  one  can  study  properties  and 
defects  in  biological  tissues  with  very  little  discomfort  to  the  subjects.  Other  potential 
applications  for  such  interrogation  techniques  are  nondestructive  damage  detection 
in  aircraft  and  spacecraft  where  very  high  frequency  electromagnetic  pulses  can  be 
used  to  detect  the  location  and  width  of  any  cracks  that  may  be  present  [BGW03]. 
Additional  applications  are  found  in  mine,  ordinance  and  camouflage  surveillance, 
and  subsurface  and  atmospheric  environmental  modeling  [BBL00]. 

Another  important  technique  employed  in  noninvasive  interrogation  of  objects  is 
the  use  of  ultrasonic  waves  in  industrial  and  medical  applications.  The  interaction  of 
electromagnetic  and  sound  waves  in  matter  has  been  widely  studied  in  acoustooptics 
[And67,  Kor97].  Brillouin  first  showed  that  electromagnetic  and  sound  waves  can 
interact  in  a  medium  and  influence  each  others  propagation  dynamics  [Bri60].  The 
response  of  the  atomic  electrons  in  a  material  medium  to  the  applied  electric  field 
results  in  a  polarization  of  the  material  with  a  resulting  index  of  refraction  that  is  a 
function  of  the  density  in  the  material.  Consequently,  the  material  density  fluctua¬ 
tions  produced  by  a  sound  wave  induce  perturbations  in  the  index  of  refraction.  Thus 
an  electromagnetic  wave  transmitted  through  the  medium  will  be  modulated  by  the 
sound  wave,  and  scattering  and  reflection  will  occur.  Conversely,  a  spatial  variation 
of  the  electromagnetic  field  intensity  and  the  corresponding  electromagnetic  stress 
lead  to  a  volume  force  distribution  in  the  medium.  This  is  called  clectrostriction,  and 
can  lead  to  sound  generation.  This  mutual  interaction  between  light  and  sound  may 
also  lead  to  instabilities  and  wave  amplification  [MI68]. 

The  work  in  this  paper  is  based  on  ideas  underlying  the  techniques  formulated  in 
[BBL00]  and  [ABR02],  In  [BBL00],  the  authors  present  an  electromagnetic  interro¬ 
gation  technique  for  general  dispersive  media  backed  by  either  a  perfect  conducting 
wall  to  reflect  the  electromagnetic  waves,  or  by  using  standing  acoustic  waves  as 
virtual  reflectors.  In  [ABR02],  this  work  was  extended  to  using  travelling  waves  as 
acoustic  reflectors.  In  both  cases  the  techniques  were  demonstrated  in  one  dimension 
by  assuming  plane  electromagnetic  waves  impinging  at  normal  incidence  on  slabs  of 
dispersive  media.  In  [BB03]  a  multidimensional  version  of  the  electromagnetic  in- 
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terrogation  problem  is  presented  for  dispersive  media  backed  by  perfect  conductor 
walls  and  the  FDTD  method  is  utilized  [Taf95,  Taf98].  ffowever,  as  we  have  noted 
above,  the  object  of  interrogation  may  not  be  backed  by  perfect  conductor  walls  in 
many  important  applications,  and  hence  one  needs  to  have  another  type  of  accessible 
interface  which  will  reflect  the  propagating  electromagnetic  pulses.  It  is  therefore 
of  great  importance  to  ascertain  (computationally  and  experimentally)  whether  the 
interrogation  ideas  of  [ABR02]  using  acoustic  pressure  waves  as  virtual  reflectors  can 
be  effective  in  multidimensional  settings.  The  goal  of  this  paper  is  the  investigation 
of  computational  and  statistical  aspects  of  this  methodology.  At  present  we  are  also 
trying  to  validate  these  proposed  techniques  by  testing  the  ideas  in  this  paper  in  an 
experimental  setup.  Our  group  has  built  an  antenna  that  we  are  using  to  study  the 
acoustooptical  interaction  in  agar-agar  in  which  the  sound  waves  are  produced  by 
means  of  a  transducer  [ABK04],  In  such  a  setup,  as  well  as  in  many  applications  one 
cannot  assume  normally  incident  plane  waves  impinging  on  the  agar,  and  hence  one 
needs  to  consider  the  more  general  case  of  oblique  incidence.  This  motivates  the  study 
of  an  electromagnetic  interrogation  technique  using  the  acoustooptical  interaction  to 
study  dielectic  properties  of  dispersive  media  in  multidimensions. 

We  begin  with  Maxwell’s  equations  for  a  first  order  dispersive  (Debye)  media 
with  Ohmic  conductivity  and  incorporate  a  model  that  describes  the  electromag¬ 
netic/acoustic  interaction.  Our  model  features  modulation  of  the  material  polariza¬ 
tion  by  the  pressure  wave  and  thus  the  behaviour  of  the  electromagnetic  pulse.  This 
approach  is  based  on  ideas  found  in  [MI68].  We  incorporate  a  pressure  dependent 
Debye  model  for  orientational  polarization  into  Maxwell’s  equations.  We  model  the 
propagation  of  a  non-harmonic  pulse  from  a  finite  antenna  source  in  free  space  im¬ 
pinging  obliquely  on  a  planar  interface  into  the  Debye  medium.  We  use  the  finite 
difference  time  domain  method  (FDTD)  [Taf95]  to  discretize  Maxwell’s  equations 
and  to  compute  the  components  of  the  electric  field  in  the  case  where  the  signal  and 
the  dielectric  parameters  are  independent  of  the  y  variable.  Figure  1  depicts  the 
antenna  and  the  Debye  dielectric  slab  geometry  that  we  will  model  in  our  problem 
with  the  infinite  dielectric  slab  perpendicular  to  the  z-axis  and  uniform  in  the  y- 
direction  and  finite  in  the  x-direction  lying  in  the  region  —  oo  <y<  oo,Xi  <  x  <  X2- 
An  alternating  current  along  the  x-direction  of  the  antenna  located  at  z  =  zc  then 
produces  an  electric  field  that  is  uniform  in  y  with  nontrivial  components  Ex  and 
E~  depending  on  ( t,x,z ).  When  propagated  in  the  xz- plane,  this  results  in  oblique 
incident  waves  on  the  dielectric  surface  in  the  xy-plane  at  z  =  Z\.  (In  [BF95]  the 
authors  use  Fourier-series  in  the  frequency  domain  to  compute  the  propagation  of 
a  time  harmonic  pulse  train  of  plane  waves  that  enter  a  dielectric  across  a  planar 
boundary  at  oblique  incidence.  They  showed  that  precursors  are  excited  by  short 
rise  time  pulses  at  oblique  incidence.  However,  the  use  of  Fourier  series  restricts 
this  approach  to  harmonic  pulses.)  The  oblique  waves  impinge  on  the  planar  inter¬ 
face  that  separates  air  and  the  Debye  medium  at  z  —  Z\  where  they  undergo  partial 
transmission  and  partial  reflection.  A  reflected  wave  travels  back  to  the  antenna  at 
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z  =  zc  and  is  recorded.  The  transmitted  part  of  the  wave  travels  through  the  Debye 
medium  and  intercepts  an  acoustic  pressure  wave  that  is  generated  in  the  dielectric 
medium.  Here  the  microwave  pulse  undergoes  partial  transmission  and  reflection. 
The  reflected  wave  travels  back  through  the  Debye  medium,  undergoes  a  secondary 
reflection/transmission  at  the  air-Debye  interface  and  is  then  recorded  at  the  antenna. 
Thus  the  antenna  records  three  pieces  of  information;  the  input  source,  the  reflection 
of  the  electromagnetic  pulse  from  the  air-Debye  interface  and  the  reflection  from  the 
pressure  wave  region.  This  data  is  recorded  at  the  center  of  the  antenna  and  will  be 
used  in  an  inverse  problem  to  estimate  the  dielectric  properties  of  the  Debye  medium. 
We  note  that  the  acoustic  wave  speed  is  many  orders  of  magnitude  smaller  than  the 
speed  of  light.  Thus  the  acoustic  pressure  wave  is  almost  stationary  relative  to  the 
microwave  pulse.  In  this  report,  we  have  assumed  that  the  acoustic  pressure  wave 
is  specified  a  priori  and  hence  we  assume  that  the  pressure  wave  is  not  modified  by 
the  electromagnetic  wave.  This  leads  to  a  much  simplified  electromagnetic/acoustic 
interaction  model.  The  dynamic  interaction  between  the  two  and  the  appropriateness 
of  this  simplifying  assumption  are  currently  being  investigated  and  will  be  the  topic 
of  a  future  report. 

An  important  computational  issue  in  modeling  two  or  three  dimensional  wave 
propagation  is  the  construction  of  a  finite  computational  domain  and  appropriate 
boundary  conditions  for  terminating  this  domain  in  order  to  numerically  simulate 
an  essentially  unbounded  domain  problem.  We  surround  the  domain  of  interest  by 
perfectly  matched  absorbing  layers  (PML’s),  thus  producing  a  finite  computational 
domain.  The  original  PML  formulation  is  due  to  J.  P.  Berenger  [Ber94]  and  was  based 
on  a  split  form  of  Maxwell’s  equations  in  which  all  the  fields  themselves  were  split 
into  orthogonal  subcomponents.  This  formulation  was  appropriately  called  the  split 
field  form.  Many  different  variations  of  the  split  field  PML  are  now  available  in  the 
literature.  Most  of  these  different  formulations,  though  equivalent  to  the  split  field 
PML,  do  not  involve  the  splitting  of  the  fields  or  splitting  of  Maxwell’s  equations. 
The  PML  formulation  that  we  use  in  this  report  is  based  on  the  anisotropic  (uniaxial) 
formulation  that  was  first  proposed  by  Sacks  et.  al.  [SKLL95]. 

We  first  present  numerical  results  for  the  forward  problem  for  two  different  test 
cases.  The  relaxation  times  of  the  Debye  media  in  these  two  test  cases  differ  by 
many  orders  of  magnitude.  We  perform  a  sensitivity  analysis  in  each  test  case  to 
determine  which  parameters  are  most  likely  to  be  accurately  estimated  in  an  inverse 
problem.  We  also  determine  which  of  the  pressure  dependent  parameters  are  the  most 
significant  in  each  of  the  test  cases.  As  we  shall  see,  the  significant  difference  in  the 
orders  of  magnitude  of  the  relaxation  times  for  the  two  test  cases  produces  inverse 
problems  with  very  different  estimation  properties. 

We  next  discuss  a  parameter  estimation  problem  in  which  we  estimate  model 
parameters  for  the  two  different  Debye  media  from  simulated  data.  We  formulate 
our  inverse  problem  using  the  method  of  least  squares.  We  compare  two  different 
algorithms  in  solving  the  parameter  estimation  problems;  the  Levenberg-Marquardt 
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method,  which  is  a  gradient  based  algorithm,  and  the  Nelder-Mead  method,  which  is 
based  on  simplices  and  is  gradient  free.  For  the  first  test  case,  these  two  algorithms 
produce  similar  results.  However,  in  the  second  test  case,  we  see  the  presence  of  many 
local  minima  and  the  Nelder-Mead  algorithm  converges  to  a  local  minima  instead  of 
the  global  minimum. 

We  also  consider  uncertainty  in  estimates  due  to  data  uncertainty.  We  present 
a  model  for  generating  noisy  data  that  we  use  in  our  inverse  problems  to  simulate 
experimental  settings  in  which  noise  enters  into  the  problem  in  a  natural  fashion.  We 
then  discuss  a  statistical  error  methodology  to  analyze  the  results  of  the  inverse  prob¬ 
lems  based  on  noisy  data.  This  statistical  error  analysis  is  used  to  obtain  confidence 
intervals  for  all  the  parameters  estimated  in  the  inverse  problem.  These  intervals 
indicate  the  level  of  confidence  that  can  be  associated  with  estimates  obtained  with 
our  methodology. 

2  Model  formulation:  The  forward  problem 

We  consider  Maxwell’s  equations  which  govern  the  electric  held  E  and  the  magnetic 
held  H  in  a  domain  hi  with  charge  density  p.  Thus  we  consider  the  system 


'  (i) 

<9D 

—  +  J  -  V  x  H  =  0,  in  (0 ,  T) 

at 

(ii) 

(9B 

—  +  V  x  E  =  0,  in  (0,  T) 

x  hi, 

(hi) 

V  •  D  =  p,  in  (0,  T)  x  H, 

(iv) 

V  •  B  =  0,  in  (0,  T )  x  H, 

(v) 

E  x  n  =  0,  on  (0,  T )  x  <9h2, 

,  (vi) 

E(0,x)  =  0,  H(0,x)  =  0, 

in  fl 

The  fields  D,  B  are  the  electric  and  magnetic  flux  densities  respectively.  All  the 
fields  in  (1)  are  functions  of  position  r  =  ( x,y,z )  and  time  t.  We  have  J  =  Jc  +  Js, 
where  Jc  is  a  conduction  current  density  and  J5  is  the  source  current  density.  We 
assume  only  free  space  (actually  the  antenna  in  our  example)  can  have  a  source 
current,  and  Jc  is  only  found  in  the  dielectric  material.  The  condition  (1,  v)  is  a 
perfectly  conducting  boundary  condition  on  the  computational  domain  which  we  will 
replace  with  absorbing  conditions  in  the  sequel.  With  perfect  conductor  conditions, 
the  computational  boundary  acts  like  a  hard  wall  to  impinging  electromagnetic  waves 
causing  spurious  reflections  that  would  contaminate  the  data  that  we  wish  to  use  in 
our  inverse  problems. 

Constitutive  relations  which  relate  the  electric  and  magnetic  fluxes  D,  B  and  the 
electric  current  Jc  to  the  electric  and  magnetic  fields  E.  H  are  added  to  these  equations 
to  make  the  system  fully  determined  and  to  describe  the  response  of  a  material  to  the 
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electromagnetic  fields.  In  free  space,  these  constitutive  relations  are  D  =  eoE,  and 
B  =  g0  H.  where  eo  and  /io  are  the  permittivity  and  the  permeability  of  free  space, 
respectively,  and  are  constant  [Jac99].  In  general  there  are  different  possible  forms 
for  these  constitutive  relationships.  In  a  frequency  domain  formulation  of  Maxwell’s 
equations,  these  are  usually  converted  to  linear  relationships  between  the  dependent 
and  independent  quantities  with  frequency  dependent  coefficient  parameters.  We  will 
consider  the  case  of  a  dielectric  in  which  magnetic  effects  are  negligible,  and  we  assume 
that  Ohm’s  law  governs  the  electric  conductivity.  Thus,  within  the  dielectric  medium 
we  have  constitutive  relations  that  relate  the  flux  densities  D,B  to  the  electric  and 
magnetic  fields,  respectively.  We  have 

(i)  D  =  e0E  +  P ID, 

(ii)  B  =yU0H.  (2) 

„  (hi)  Jc  =  crE ID, 

In  (2),  If)  denotes  the  indicator  function  on  the  Debye  medium.  Thus,  Jc  =  0  in 
air.  The  quantity  P  is  called  the  electric  polarization.  It  is  equal  to  zero  in  air  and 
is  nonzero  in  the  dielectric.  Electric  polarization  may  be  defined  as  the  electric  held 
induced  disturbance  of  the  charge  distribution  in  a  region.  This  polarization  may 
have  an  instantaneous  component  as  well  as  ones  that  do  not  occur  instantaneously; 
the  latter  usually  have  associated  time  constants  called  the  relaxation  times  which  are 
denoted  by  r.  We  let  the  instantaneous  component  of  the  polarization  to  be  related 
to  the  electric  held  by  means  of  a  dielectric  constant  eoX  and  denote  the  remainder 
of  the  electric  polarization  as  Pr.  We  have 


P  =  Pi  +  Pr  =  eoyE  +  Pr, 
and  hence  the  constitutive  law  (2,  i)  becomes 

D  =  eoerE  +  Pr, 


(3) 

(4) 


where  er  —  (1  +  x)  is  the  relative  permittivity  of  the  dielectric  medium. 

To  describe  the  behaviour  of  the  media’s  macroscopic  electric  polarization  Pr, 
we  employ  a  general  integral  equation  model  in  which  the  polarization  explicitly  de¬ 
pends  on  the  past  history  of  the  electric  held.  This  model  is  sufficiently  general  to 
include  microscopic  polarization  mechanisms  such  as  dipole  or  orientational  polariza- 
tion(Debye),  as  well  as  ionic  and  electronic  polarization(Lorentz),  and  other  frequency 
dependent  polarization  mechanisms  [And67]  as  well  as  multiples  of  such  mechanisms 
represented  by  a  distribution  of  the  associated  time  constants  (e.g.,  see  [BG04]).  The 
resulting  constitutive  law  can  be  given  in  terms  of  a  polarization  or  displacement 
susceptibility  kernel  g  as 


PrC t,x,z) 


g(t  —  s,  x,  z)E(s,  x,  z)ds, 


(5) 
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which,  in  the  case  (Debye)  of  interest  to  us  here,  can  be  written  as 

g(t)  =  e-br£o(^e  oo). 

r 

Such  a  Debye  model  can  be  represented  in  differential  form  as 

'tPr  +  Pr  —  Co(es  —  Coo)E) 

D  =  eofiooE  +  Pr, 


(6) 

(7) 


inside  the  dielectric,  whereas  Pr  =  0,  =  1  in  air,  and  =  er  in  the  dielectric. 

We  will  henceforth  denote  Pr  by  P.  In  equations  (5)  and  (7),  the  parameters  es 
and  Coo  denote  the  static  relative  permittivity,  and  the  value  of  permittivity  for  an 
extremely  high(«  oo)  frequency  field,  respectively.  The  quantity  is  called  the  infi¬ 
nite  frequency  permittivity.  Biological  cells  and  tissues  display  very  high  values  of  the 
dielectric  constants  at  low  frequencies,  and  these  values  decrease  in  almost  distinct 
steps  as  the  excitation  frequency  is  increased.  This  frequency  dependence  of  permit¬ 
tivity  is  called  dielectric  dispersion  and  permits  identification  and  investigation  of  a 
number  of  completely  different  underlying  mechanisms.  This  property  of  dielectric 
materials  makes  the  problem  of  parameter  identification  an  important  as  well  as  very 
useful  topic  for  investigation. 

We  will  modify  system  (1)  and  the  constitutive  laws  (2)  by  performing  a  change 
of  variables  that  renders  the  system  in  a  form  that  is  convenient  for  analysis  and 
computation.  From  (1,  i)  we  have 


d_ 

dt 


—  VxH  =  — Js,  in  (0,  T )  x  fl 
8 


(8) 


Next,  we  define  the  new  variable 


D(f,  x)  =  D(f,  x)  +  /  Jc(s,x)  ds. 


Using  definition  (9)  in  (8)  and  (1)  we  obtain  the  modified  system 


(i)  ^  -  VxH  =  —  J,  in  (0,T)  xfi, 

dB 

(ii)  — — h  V xE  =  0  in  (0,  T)  x  ff, 

(y  L 

(iii)  V  •  D  =  0  in  (0,  T )  x  fl, 

(iv)  V  •  B  =  0  in  (0,  T )  x  U, 

(v)  E  x  n  =  0  on  (0,  T)  x  <9h2, 

(vi)  E(0, x)  =  0,  H(0,x)  =  0  in  Q. 


(9) 


(10) 


We  will  henceforth  drop  the 'symbol  on  D.  We  note  that  equation  (10,  iii)  follows 
from  the  continuity  equation  ^  _|_  y  .  j  —  o,  the  assumption  that  p(0)  =  0,  and  the 
assumption  that  V  •  Js  =  0  (in  the  sense  of  distributions).  The  modified  constitutive 
laws  (2)  after  substitution  of  Ohm’s  law  and  (9)  are 

(i)  D(t,x)  =  e0eoo(x)E(f,x)  +  PR  +  J  <r(x)E(s,  x)  ds, 

,  (ii)  B  =  p0H. 

with  PR,  the  solution  of 

'tPr  +  Pr  —  ^0  (es  —  eoc)E.  (12) 

3  An  anisotropic  perfectly  matched  layer  absorb¬ 
ing  medium 

In  one  dimension  exact  absorbing  boundary  conditions  can  be  constructed  for  ter¬ 
minating  a  computational  domain.  In  two  and  three  dimensions,  however,  there  are 
no  exact  absorbing  boundary  conditions  available.  In  such  a  case  we  can  either  use 
approximate  absorbing  boundary  conditions  of  some  finite  order  (usually  one  or  two), 
or  we  can  surround  our  computational  domain  by  absorbing  layers.  In  this  report  we 
use  the  second  approach  and  construct  perfectly  matched  layers  [Ber94]  surrounding 
the  domain  of  interest.  These  layers  have  been  demonstrated  to  have  many  orders  of 
magnitude  better  absorbing  capabilities  than  the  approximate  absorbing  boundary 
conditions  [Taf98].  In  this  section  we  will  first  construct  PML’s  that  can  be  matched 
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to  free  space  and  then  we  will  explain  how  this  construction  can  be  generalized  to  the 
case  of  a  Debye  medium. 

The  source  current  JA,  is  zero  inside  the  PML  regions.  So  we  will  neglect  Js  in 
this  section.  We  model  the  PMLs  as  lossy  regions  by  adding  artificial  loss  terms  to 
Maxwell’s  equations.  This  approach  is  based  on  using  anisotropic  material  properties 
to  describe  the  absorbing  layer  [SKLL95,  Ged96,  Rap95].  We  start  with  Maxwell’s 
equations  in  the  most  general  form 


f  dD 

~di 

<9B 

~dt 


VxH  —  Je, 

VxE  Jrri) 


V-D  =  0, 
,  V-B  =  0, 


(13) 


where  Je  and  Jm  are  electric  and  magnetic  conductivities,  respectively.  Inside  the 
computational  domain  we  will  have  Je  =  Jc  and  Jm  =  0.  We  will  derive  a  PML 
model  in  the  frequency  domain  and  then  transform  the  corresponding  equations  to 
the  time  domain  by  taking  the  inverse  Fourier  transforms  of  the  frequency  domain 
equations.  To  this  end,  we  consider  the  time-harmonic  form  of  Maxwell’s  equations 
(13)  given  by 

jcuD  =  VxH  —  Je, 
jcnB  —  V  xE  Jm, 

(14) 

V-B  =  0, 

,  V-D  =  0, 


where  for  every  held  vector  V,  V  denotes  its  Fourier  transform, 
constitutive  laws 


(  B  =  [/i]H, 


We  then  have  the 


D  =  [e]E. 

J  m  [Fm]H, 
,  Je  =  We] E. 


(15) 


In  (15)  the  square  brackets  indicate  a  tensor  quantity.  Thus  we  allow  for  anisotropy  in 
our  media  by  permitting  the  material  parameters  to  be  tensors.  Note  that  when  the 
density  of  electric  and  magnetic  charge  carriers  in  the  medium  is  uniform  throughout 
space,  then  V  •  Je  =  0  and  V  •  Jm  =  0.  We  define  the  new  tensors 


-  W\  + 


e  =  e  + 


Wm\ 

.w 

We] 


(16) 
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Using  the  definitions  (16)  we  define  two  new  constitutive  laws  that  are  equivalent  to 
(15),  given  by 

Bnew  =  [/i]H, 

A  A  (17) 

„  Dnew  —  [e]E. 

Using  (17)  in  (14),  we  can  write  Maxwell’s  equations  in  time-harmonic  form  as 

'  jcdDnew  =  VxH. 
jcaBnew  =  V  xE, 

(18) 

V  •  Duew  =  0, 

,  V  •  Bnew  =  0. 

The  split-field  PML  introduced  by  Berenger  [Ber94]  is  a  hypothetical  medium 
based  on  a  mathematical  model.  In  [MP95]  Mittra  and  Pekel  showed  that  Berenger’s 
PML  was  equivalent  to  Maxwell’s  equations  with  a  diagonally  anisotropic  tensor 
appearing  in  the  constitutive  relations  for  D  and  B.  For  a  single  interface,  the 
anisotropic  medium  is  uniaxial  and  is  composed  of  both  the  electric  and  magnetic 
permittivity  tensors.  This  uniaxial  PML  (UPML)  formulation  performs  as  well  as 
the  original  split-field  PML  while  avoiding  the  nonphysical  field  splitting.  As  will  be 
shown  below,  by  properly  defining  a  general  constitutive  tensor  [S'],  we  can  use  the 
UPML  in  the  interior  working  volume  as  well  as  the  absorbing  layer.  This  tensor 
provides  a  lossless  isotropic  medium  in  the  primary  computation  zone,  and  individual 
LTPML  absorbers  adjacent  to  the  outer  lattice  boundary  planes  for  mitigation  of 
spurious  wave  reflections.  The  fields  excited  within  the  LTPML  are  also  plane  wave  in 
nature  and  satisfy  Maxwell’s  curl  equations.  The  derivation  of  the  PML  properties 
of  the  tensor  constitutive  laws  is  also  done  directly  by  Sacks  et  al.  in  [SKLL95]  and 
by  Gedney  in  [Ged96].  We  follow  the  derivation  by  Sacks  et  al.  here.  We  begin 
by  considering  planar  electromagnetic  waves  in  free  space  incident  upon  a  PML  half 
space.  Starting  with  the  impedance  matching  assumption,  i.e.,  the  impedance  of  the 
layer  must  match  that  of  free  space:  eg  Vo  =  [e]_1V]  we  have 

—  =  —  =  [S]  =  diag{ai,  ai,  a3}.  (19) 

e0 

Hence,  the  constitutive  parameters  inside  the  PML  layer  are  [e]  =  eo[S]  and  [p]  = 
yUoV],  where  [S']  is  a  diagonal  tensor.  We  next  consider  plane  wave  solutions  of  the 
form 

V(x,  t)  =  V(x)  exp(j(u;f  —  k  •  x)),  (20) 

for  all  field  vectors  V,  to  the  time-harmonic  Maxwell’s  equations  with  the  diagonally 
anisotropic  tensor,  where  k  =  ( kx ,  ky ,  kz )  is  the  wave  vector  of  the  planar  electromag¬ 
netic  wave  and  x  =  (x,y,z).  The  dispersion  relation  for  waves  in  the  PML  are  found 
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to  be 


(21) 


kl 


+ 


K 


+ 


kl 


—  /Cq  =  Ci/^oU)  =  — 2~. 


cl  2CI3  aia3  Q1O2 

where,  c0  is  the  speed  of  light  in  free  space. 

To  further  explain  our  UPML  formulation,  we  may,  without  loss  of  generality, 
consider  a  PML  layer  which  fills  the  positive  x  half-space  and  plane  waves  with  wave 
vectors  in  the  xz-  plane  (ky  =  0).  Let  9,  be  the  angle  of  incidence  of  the  plane  wave 
measured  from  the  normal  to  the  surface  x  —  0.  The  standard  phase  and  magnitude 
matching  arguments  at  the  interface  yield  the  following  generalization  of  Snell’s  law 


y/ai03sm9t  =  sin  9i: 


(22) 


where  9t  is  the  angle  of  the  transmitted  plane  wave.  Matching  the  magnitudes  of  the 
electric  and  magnetic  fields  at  the  interface,  x  =  0,  we  obtain  the  following  values  of 
the  reflection  coefficients  for  the  TE  and  the  TM  modes 

(  a  /03  Q 

cos  t h~\  —  cos  “t 

V  a2 


R 


TE 


R 


TM 


cos  9{  +  \  —  cos  9t 

V  «2 

[&■■ 3  n  a 

—  COS  9t  —  COS  9 1 

cos  9{  +  \  —  cos  9t 

V  a2 


(23) 


From  (23)  we  can  see,  that  by  choosing  a3  =  02  =  a  and  ^/cRa^  =  1,  the  interface 
is  completely  reflectionlcss  for  any  frequency,  angle  of  incidence,  and  polarization. 
Using  (17)  and  (19),  we  thus  find  the  constitutive  laws  for  the  perfectly  matched 
layer  are 


Dr 

B, 


eo[S}E, 

ho[^]H. 


where  the  tensor  [S']  is 


[5]  = 


x  1 
0 
0 


0 

a 

0 


0 

0 

a 


(24) 


(25) 


The  perfectly  matched  layer  is  therefore  characterized  by  the  single  complex  number 
a.  Taking  it  to  be  the  constant  a  =  7  —  i/3,  and  substituting  into  the  dispersion 
relation  (21),  we  obtain  the  following  expression  for  the  electric  held  inside  of  the 
PML 

E (x,y,z)  =  E0  exp(— ko/3  cos9tx)  exp(—  j/c0(7  cos  9tx  +  sin  6^2))  exp(jcuf).  (26) 

Hence  we  can  see  that  7  determines  the  wavelength  of  the  wave  in  the  PML  and  for 
(3  >  0,  the  wave  is  attenuated  according  to  the  distance  of  travel  in  the  x  direction. 
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Remark  1  To  match  the  anisotropic  PML  along  a  planar  boundary  to  a  dispersive 
isotropic  half-space,  we  will  need  to  modify  the  above  construction  for  the  phasor- 
domain  constitutive  relation 

D  =  €o6r(u)it,  (27) 

where  er(oj)  is  the  frequency  dependent  relative  permittivity  of  the  Debye  medium, 
which  will  be  defined  in  a  standard  manner  in  (40)  below. 

3.1  Implementation  of  the  uniaxial  PML 

To  apply  the  perfectly  matched  layer  to  electromagnetic  computations,  we  replace  the 
half  infinite  layer  by  a  layer  of  finite  depth  backed  with  a  more  conventional  boundary 
condition,  such  as  a  perfect  electric  conductor  (PEC).  This  truncation  of  the  layer  will 
lead  to  reflections  generated  at  the  PEC  surface,  which  can  propagate  back  through 
the  layer  to  re-enter  the  computational  region.  In  this  case,  the  reflection  coefficient 
R ,  is  now  a  function  of  the  angle  of  incidence  6,  the  depth  of  the  PML  6,  as  well  as 
the  parameter  a  in  (25).  Thus,  this  parameter  a  for  the  PML  is  chosen  in  order  for 
the  attenuation  of  waves  in  the  PML  to  be  sufficient  so  that  the  waves  striking  the 
PEC  surface  are  negligible  in  magnitude.  Perfectly  matched  layers  are  then  placed 
near  each  edge  (or  face  in  the  3D  case)  of  the  computational  domain  where  a  non¬ 
reflecting  condition  is  desired.  This  leads  to  overlapping  PML  regions  in  the  corners 
of  the  domain.  As  shown  in  [SKLL95],  the  correct  form  of  the  tensor  which  appears 
in  the  constitutive  laws  for  these  regions  is  the  product 

[S]  =  [S],[S]„[S]„  (28) 

where  component  [S] a  in  the  product  (28)  is  responsible  for  attenuation  in  the  a 
direction,  for  a  =  x,y,z  (see  Figure  2).  All  three  of  the  component  tensors  in  (28) 
are  diagonal  and  have  the  forms 

s^1  0  0  sy  0  0  sz  0  0 

[S\x=  0  sx  0  ;[S}y=  0  0  ;[£],=  0  sz  0  .  (29) 

0  0  sx  0  0  sy  0  0  S71 

In  the  above  sx,sy,sz  are  analogous  to  the  complex  valued  parameter  a  encountered 
in  the  Section  3  analysis  of  the  single  PML  layer.  Thus,  sa  governs  the  attenuation  of 
the  electromagnetic  waves  in  the  a  direction  for  a  =  x,y,  z.  When  designing  PML’s 
for  implementation,  it  is  important  to  choose  the  parameters  sQ  so  that  the  resulting 
frequency  domain  equations  can  be  easily  converted  back  into  the  time  domain.  The 
simplest  of  these  [Ged96]  which  we  employ  here,  is 

sa  —  1  4 - — ,  where  aa  >  0,  a  =  x,  y,  z.  (30) 

jwe0 
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Figure  2:  PML  layers  surrounding  the  domain  of  interest.  In  the  corner  regions  of  the 
PML,  both  ax  and  az  are  positive  and  the  tensor  [S']  is  the  product  [S']x[5,]z.  In  the 
remaining  regions  only  one  of  oz  (left  and  right  PML’s)  or  crx  (top  and  bottom  PML’s) 
are  nonzero  and  positive.  The  tensor  [S'],  is  thus  either  [S']*  or  [S']2  respectively.  The 
PML  is  truncated  by  a  perfect  electric  conductor  (PEC). 


The  PML  interface  represents  a  discontinuity  in  the  conductivities  cra.  To  reduce 
the  numerical  reflections  caused  by  these  discontinuous  conductivities,  the  aa  are 
chosen  to  be  functions  of  the  variable  a  (e.g.,  ax  is  taken  to  be  a  function  of  x  in 
the  [S]x  component  of  the  PML  tensor).  A  choice  of  these  functions  so  that  aa  =  0, 
i.e.,  sa  —  1  at  the  interface  yields  the  PML  a  continuous  extension  of  the  medium 
being  matched  and  reduces  numerical  reflections  at  the  interface.  If  one  increases 
the  value  of  cra  with  depth  in  the  layer,  one  obtains  greater  overall  attenuation  while 
minimizing  the  numerical  reflections.  Gedney  [Ged96]  suggests  a  conductivity  profile 


aa(a) 


8m 


a  =  x,y,  z. 


(31) 


where  5  is  the  depth  of  the  layer,  a  =  ao  is  the  interface  between  the  PML  and 
the  computational  domain,  and  m  is  the  order  of  the  polynomial  variation.  Gedney 
remarks  that  values  of  m  between  3  and  4  are  believed  to  be  optimal.  For  the 
conductivity  profile  (31),  the  PML  parameters  can  be  determined  for  given  values  of 
m,  5,  and  the  desired  reflection  coefficient  at  normal  incidence  Rq,  as 


^"max 


(m  +  1)  ln(l/i?0) 
2  r/i8 


(32) 


r]j  being  the  characteristic  wave  impedance  of  the  PML.  Emperical  testing  suggests 
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that,  for  a  broad  range  of  problems,  an  optimal  value  of  crmax  is  given  by 


m  +  1 

a°pt  “  150ttAox/£;’ 


(33) 


where  Aa  is  the  space  increment  in  the  a  direction  and  er  is  the  relative  permittivity 
of  the  material  being  modeled.  In  the  case  of  free  space  er  —  1.  Since  we  desire  to 
match  the  PML  to  both  free  space  as  well  as  the  Debye  medium,  we  will  use  er  to  be 
the  average  value 

+  too)  .  (34) 


3.2  Reduction  to  two  dimensions 


From  the  time-harmonic  Maxwell’s  curl  equations  in  the  PML  (18)  and  (24),  Ampere’s 
and  Faraday’s  laws  can  be  written  in  the  most  general  form  as 

r  j^0[s]H  —  — X  xE,  (Faraday  s  Law) 


^  jca[S]D  =  VxH  —  Js,  (Ampere’s  Law) 


In  (35),  [5]  is  the  diagonal  tensor  defined  via  (28),  (29)  and  (30).  In  the  presence  of 
this  diagonal  tensor,  a  plane  wave  is  purely  transmitted  into  the  uniaxial  medium. 
The  tensor  [S']  is  no  longer  uniaxial  by  strict  definition,  but  rather  is  anisotropic. 
However,  the  anisotropic  PML  is  still  referenced  as  uniaxial  since  it  is  uniaxial  in 
the  nonoverlapping  PML  regions.  Let  Q  =  X  x  Z  denote  the  computational  domain 
along  with  the  absorbing  layers.  We  partition  the  interval  X  into  disjoint  closed 
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intervals  as  X  =  XBpmi  U  X0  U  Xrpmi  and  the  interval  Z  into  disjoint  closed  intervals 
as  ^Lpmi  uzau  Zjy  U  ^Rpmi ,  as  seen  in  Figure  3.  The  computational  domain  of  interest 
is  the  region  X0  x  {Za  U  ZDj,  where  X0  x  Za  denotes  air  and  X0  x  ZD  denotes  the 
Debye  medium.  The  buffer  region  T>  \  (X0  x  {Za  U  ZD  })  contains  the  absorbing 
layers  (PML’s).  The  PML’s  are  backed  by  a  perfect  conductor  where  the  boundary 
condition  n  x  E  =  0  is  used.  We  construct  perfectly  matched  layers  in  regions 
where  (x,z)  G  (ATspmi  UlTpmi)  x  (^Lpmi  U  Zptpmi).  To  obtain  a  two-dimensional 
model,  we  make  the  assumption  that  the  signal  and  the  dielectric  parameters  are 
independent  of  the  y  variable.  Also,  an  alternating  current  along  the  x  direction,  as 
in  (54),  then  produces  an  electric  field  that  is  uniform  in  y  with  nontrivial  components 
Ex  and  Ez  depending  on  (t,  x,  z)  which,  when  propagated  in  the  xz  plane  result  in 
oblique  incident  waves  on  the  dielectric  surface  in  the  xy  plane  at  z  =  zc.  With  these 
assumptions,  Maxwell’s  equations  reduces  to  a  two  dimensional  system  for  solutions 
involving  the  Ex  and  Ez  components  for  the  electric  field  and  the  component  Hy  for 
the  magnetic  field,  which  is  refered  to  as  the  TEy  mode.  In  this  case,  we  have  <ry  =  0 
and  sy  =  1  in  the  UPML.  From  (11,  i)  and  (7)  expressed  in  the  frequency  domain  we 
have  the  constitutive  relation 


Dl  .  eoo  .  ® 

-  to  Poo  +  Z - - b  7 - 

l+jcur  jue0 


E. 


Rescaling  the  electric  field  as 


~E=\-  E. 
V  Mo 


(36) 


(37) 


we  can  write  the  time-harmonic  frequency  domain  Maxwell’s  equations  (35)  in  the 
uniaxial  medium  as 


(Tz(z) 
j^o 

o~x(x ) 

j^e0 

/...x  ■  hl  0x(x) 

(ill)  jw  1+  - - 

J^eo 


(i)  j u  ^1  + 

(ii)  ]uj(l  + 


1  + 


1  + 


-l 


<*x{x) 
j^e0 

<d(A>  V*  A  „  dHv 


A  ,dHy  T  .. 
Dx  =  -c0(- +  Js  ■  x), 


J^e0 


Ez  Cq 


1  ,  a*(z)  \  u 

1  +  V -  Hy  —  Co 

J^e0 


dx 

BE,  dEr 


(38) 


dx  dz 


In  the  above  Co  =  3.0  x  108  m/s  is  the  speed  of  light  in  vacuum  and  D  =  (Dx,  DZ)T 
is  defined  as 

D  =  e*  •  E,  (39) 

with 

*  i  ^ S  e00  .  ®  (  AC\\ 

er  ~  eoo  +  Z — TV - b  7 - •  (40) 

1+JCCT  JCCCo 

To  avoid  a  computationally  intensive  implementation,  we  define  suitable  constitutive 
relationships  that  facilitate  the  decoupling  of  the  frequency  dependent  terms  [Taf98]. 
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To  this  end,  we  introduce  the  fields 


D*  =  s~LD 
n*  =  s~lD 

^  z  az 


H*y  =  SZHy. 


(41) 


Again,  dropping  the 'symbol  on  D  and  transforming  the  frequency  domain  equations 
to  the  time  domain,  we  obtain  the  system 

„  ax(x)  rr.  / dE z  dEx 

d,H ,  +  — =  d,H‘ 


11 


Co 


(ii.)  ajrm  +  °-f-irm  =  -e^ 

ov)  + 

(V)  dtD,  +  rcoM» 


+  Js  •  x), 


eo 

(vi)  dtDz  =  dtD*z  + 


<Jz(z) 


Co 


<9x 


(42) 


with 

D(t)  =  e0OE(t)  +  P(t)  +  C(t),  (43) 

where  D  =  (DX1DZ)T,  and  inside  of  the  dielectric  the  polarization  P  =  (Px,  PZ)T 
satisfies 

rP  +  P  =  (es  —  e00)E,  (44) 

while  the  conductivity  term  C  =  (Cx,Cy)T  satisfies 

C  =  (cr/e0)E  .  (45) 


inside  of  the  dielectric.  Outside  the  dielectric  we  have  P  =  C  =  0.  We  will  assume 
zero  initial  conditions  for  all  the  fields. 


3.3  Pressure  dependence  of  polarization 

We  consider  a  pressure  dependent  Debye  model  for  orientational  polarization  first 
developed  in  [ABR02],  We  assume  that  the  material  dependent  parameters  in  the 
differential  equation  (44)  for  the  Debye  medium  depend  on  pressure  p,  i.e., 

t(p)P  +  P  =  (ea(p)  -  eoo(p))E.  (46) 
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We  suppose  as  a  first  approximation  that  each  of  the  pressure  dependent  parameters 
can  be  prepresented  as  a  mean  value  plus  a  perturbation  that  is  proportional  to  the 
pressure 

(i)  r(p)  =  t0  +  nTp, 

(ii)  es(p)  =  esfi  +  Ksp,  (47) 

(hi)  Coo  (p)  =  600,0  +  hiocP- 

This  can  be  considered  as  a  truncation  after  the  first-order  terms  of  power  series 
expansion  of  these  functions  of  p.  Thus  (44)  is  written  as 

(t0  +  Krp) P  +  P  =  ((es,o  -  £00,0)  +  (ks  -  Koo)p)E.  (48) 

This  model  features  modulation  of  the  material  polarization  by  the  pressure  wave  and 
thus  the  behaviour  of  the  electromagnetic  pulse.  We  will  assume  that  the  acoustic 
pressure  wave  is  specified  a  priori.  The  pressure  p  will  be  assumed  to  have  the  form 

pit,  z )  =  /(Z2,22+ Xp)  sin  (up[t  +  ,  (49) 

where  Z2  €  Z d  with  z\  <  z2-  The  terms  up,  \p  and  cp  denote  the  acoustic  frequency, 
wavelength  and  speed  respectively.  We  note  that  outside  the  pressure  region,  the 
pressure  coefficients  kt,  ks  and  are  all  zero  and 


7”  7*0)  ^ s  ^s,  0?  ^00  ^00, 0* 

(50) 

From  (43) 

c 

4)  Dx(t)  =  (Coo^  +  ^oo  p)Ex(t)  +  Px(t)  +  Cx(t), 

,ii)  Dz{t)  =  (eoo,o  +  K00  p)Ez(t)  +  Pz{t)  +  Cz(t), 

(51) 

with 

( !«> 

(To  T  KrP)Px  T  Px  ii^sfl  £00, 0)  T  {ks  Koo)p)Ex, 

(to  T  KtP)Pz  T  Pz  ((£«,o  £00, 0)  T  i^s  ^oo)p) P zi 

(52) 

and 

(  (i)  Cx  =  ( cf/e0)Ex , 

[  (ii)  Cz  =  (a/eo)Ez. 

(53) 

3.4 

Source  term 

The  source  term  JA.  will  model  an  infinite  (in  the  y  direction)  antenna  strip,  finite  in 
the  x  direction  (between  X\  and  x2),  lying  in  the  z  =  zc  plane  in  free  space,  and  we 
assume  the  signal  is  polarized  with  oscillations  in  the  xz  plane  only,  with  uniformity 
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in  the  y  direction  (see  again  Figure  1).  We  therefore  take 


J s(t,x,z)  =  I(X1jX2)(x)6(z  -  zc)  sm(uc(t  -  3f0))exp 


t  —  3f0 

1 

to 

2 

to  i  ^ c 

UJd 


2i t/c  rad/sec. 


(54) 


where  I(Xi,x2)  is  the  indicator  function  on  the  interval  (xi,X2)  G  X0,  S(z  —  zc)  is 
the  Dirac  measure  centered  at  z  =  zCl  and  x  is  the  unit  vector  in  the  x  direction. 
The  Fourier  spectrum  of  this  pulse  has  even  symmetry  about  the  frequency  fc.  The 
pulse  is  centered  at  time  step  3 10  and  has  a  1/e  characteristic  decay  of  t0  time-steps 
[Taf95] .  We  will  specify  the  values  of  all  the  different  quantities  involved  in  (54) 
when  we  define  the  particular  problems  that  we  will  solve.  The  source  (54)  radiates 
numerical  waves  having  time  waveforms  corresponding  to  the  source  function.  In 
the  discretized  model,  described  below  the  radiated  numerical  wave  will  propagate  to 
the  Debye  medium  and  undergo  partial  transmission  and  partial  reflection,  and  time 
stepping  can  be  continued  until  all  transients  decay. 


3.5  Space  and  time  discretization 

We  will  use  the  FDTD  algorithm  [Taf95]  to  discretize  Maxwell’s  equations.  The 
FDTD  algorithm  uses  forward  differences  to  approximate  the  time  and  spatial  deriva¬ 
tives.  The  FDTD  technique  rigorously  enforces  the  vector  held  boundary  conditions 
at  interfaces  of  dissimilar  media  at  the  time  scale  of  a  small  fraction  of  the  impinging 
pulse  width  or  carrier  period.  This  approach  is  very  general  and  permits  accurate 
modeling  of  a  broad  variety  of  materials  ranging  from  living  human  tissue  to  radar 
absorbers  to  optical  glass.  The  system  of  equations  (42)  can  be  discretized  on  the  stan¬ 
dard  Yee  lattice.  The  domain  D  is  divided  into  square  cells  of  length  h  =  Ax  =  A z, 
h  being  the  spatial  increment.  The  degrees  of  freedom  on  an  element  are  shown  in 
Figure  4.  The  electric  held  degrees  of  freedom  are  at  the  midpoints  of  edges  of  the 
squares  and  the  magnetic  degrees  of  freedom  are  at  the  centers  of  cells.  The  time 
interval  over  which  time  stepping  is  done  is  divided  into  subintervals  of  equal  size 
using  the  time  increment  At.  We  perform  normal  leapfrogging  in  time  and  the  loss 
terms  are  averaged  in  time.  This  leads  to  an  explicit  time  stepping  scheme.  In  this 
case,  a  stability  condition  (CFL)  has  to  be  satished  [Taf95]  in  order  to  obtain  a  well 
posed  computational  scheme.  The  time  step  At  and  the  spactial  increment  h  for  the 
FD-TD  scheme  in  non-dispersive  dielectics  are  related  by  the  condition 


hex 


Co  At  1 

h  v/2 


(55) 


The  number  t/cn  is  called  the  Courant  number.  We  hx  the  value  of  ?7cn  =  1/2. 
In  [Pet94],  the  author  established  that  several  extensions  of  the  FD-TD  for  Debye 
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Edge  2 
Ez 


Figure  4:  Degrees  of  freedom  of  E  and  H  on  a  rectangle. 


media  satisfy  the  stability  restriction  of  the  standard  FD-TD  scheme  in  nondispersive 
dielectrics.  He  has  also  shown  that  the  extensions  do  not  preserve  the  non-dissipative 
character  of  the  standard  FD-TD  scheme.  The  extended  difference  systems  are  more 
dispersive  than  the  standard  FD-TD,  and  their  accuracy  depends  strongly  on  how 
well  the  chosen  timestep  resolves  the  shortest  timescale  in  the  problem  regardless  of 
whether  it  is  the  incident  pulse  timescale,  the  medium  relaxation  timescale,  or  the 
medium  resonance  timescale. 

For  any  field  variable  V,  we  use  the  notation  Vfm  =  V (rAt,  Ih,  mh ),  where  r,  /,  m  G 
N  or  r  =  r0  ±  1/2,  l  —  l0  ±  1/2,  m  —  m0  ±1/2,  with  r0,  Z0,  m0  G  N.  With  the  degrees  of 
freedom  as  shown  in  Figure  4,  the  Maxwell  curl  equations  are  discretized  as  follows. 
Given  the  values  of  the  electric  field  at  time  riAt  and  the  values  of  the  magnetic  field 
at  time  (n  —  1/2) At,  we  then  evaluate  the  magnetic  field  at  time  (n+  1/2) At  and  the 
electric  held  at  time  (n  +  l)At.  We  have  for  H*  =  H* 


20 


Solving  (56)  for  H* \™+i/2k+i/2  we  h&ve  ^ie  update  equation 

TT*,n+i/2  _  f 2ep  —  <Jx(i  +  l/2)At\  n_i/2  ( _ 2e0CoAt _ 

li+i/2,fc+i/2  -  ^260  +  ^^ +  1/2) At )  l*+l/2,fe+l/2  \^2e0  +  ax(i  +  l/2)At 


(E,\ 

\n  Z71 

li+l,AH-l/2  £L'Z\ 

71  /T* 

li,fc+l/2  ^1 

In  ZT1 

2+1/2, fc+1  ^1 

71  \ 

i+1/2, k  \ 

\ 

Arc 

Az 

(57) 


Similarly  we  discretize  the  update  equations  for  the  other  held  variables  in  (42).  The 
update  equation  for  Hy  =  H  is 


TT  I  71+1/2 

n  U+l/2,fc+l/2 

+ 


f  2ep  —  crz(k  +  l/2)At\  n_!/2 

\2e0  +  crz{k  +  l/2)At  y  |i+i/2,fc+i/2 


2en 


2eo  +  az{k  +  l/2)At 


x 


(  LX*  1 71+1/2  rr*  in— 1/2 

^  U+l/2,/c+l/2  _  U+l/2,fc+l/2 


The  update  equation  for  D*x  is 


1Jx\i+l/2,k  ~ 


f  2e0c0At  \  ^ 
\2e0  +  az(k)At)  * 


{2e0  -  az(k)At\ 

\2e0  +  az(k)At )  A+1/2’fe 


77 


I  71+ 1/2 


27  U+l/2,fc+l/2 


H, 


177+1/2 


27  U+l/2,fc— 1/2 


A* 


+  «/< 


s, 2+1/2, /c 


5 


(58) 


(59) 


where  the  electromagnetic  input  source  is  chosen  to  have  the  form  J++1/,2  k  =  J++1 /2  fcx- 
The  update  equation  for  Dx  is 


D 


171+1  _  p<  yn  I 

x\i+l/2,k  lyx\i+l/2,k' 


The  update  equation  for  D*  is 


2e0  +  ax{i  +  1/2)A+ 


2eo 

2e0  -  ax(i  +  l/2)At\ 


)D 


2en 


x\i+l/2,k 

*  |7i 

x \i-\-l/2,k  • 


7~j*  in+1 
iyz\i,k+l/2 


2e0  -  ax{i)At 
2e0  +  crx(i)At 

2e0c0At 
2e0  +  ax(i)At 

Finally  the  update  equation  for  Dz  is 


7-)*  |  n 

-Uz\i,k+l/2 


H, 


x 


171+1/2 

27U+l/2,fc+l/2 


—  H, 


171+1/2 

yli-i/2,fc+i/2 


Arc 


D 


|  n+1 

^U,fc+l/2 


=  L> 


2  ll,fc+l/2 


+ 


2eo  +  &z(k  +  l/2)At\ 


2eo  / 

2e0  -  crz{k  +  l/2)At\ 


r^*  |n+l 
U,AH-l/2 


2en 


1  n*in 
I  1 /2 * 


(60) 


(61) 


(62) 
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For  the  update  of  Ex  and  Ez  we  will  combine  equations  (51)-(53).  From  (53)  we  have 
the  update  for  Cx  and  Cz  as 


'  r  |n+l  _  r  in  ,  ( 

^  X  I 


cz\n+1  =  cz\n  + 

From  (52)  we  have  the  update  for  Px  as 

-l 


2eo 

aAt\  (P|”+1  +  PI") . 


(63) 


2e 


px  r+i  =  1 1+ 

H1+ 

+  (i  + 


At 


1  - 


At 


Pr 


2  (r0  +  KrPn+l )  J  V  2(r0  +  nTp 1 

A t  \  ~ 1  /  At((es,0  -  600,0)  +  (ks  -  Koo)pn) 


2(r0  +  KTpn+l)J  V  2(ro  +  «Tpn) 

At  \  -1  /  At((ea,0  -  600,0)  +  (/ca  -  Koo)Pn+1)\  |w+1 


2(r0  +  6cTpn+1)  /  V  2(r0  +  nTpn+l) 

Similarly,  from  (52)  we  have  the  update  for  Pz  as 

At  \_1  ( .  At 


Er\n 
Ex 


(64) 


P  "+1  =11  + 

+  2(tq  +  KTpn+1) 


1  - 


2(r0  +  kt  pr 


P, 


H1+ 


At 


2  (r0  +  KTpn+1 ) 
At 


-1 


At ( ( 6.s,o  6oo,o)  "F  (k.s  ^oo)V  ) 

2  (Tq  +  Krpn) 


Ex 


(65) 


2(r0  +  KTpn+1' 
From  (51)  we  have 


1  ( At((es,0  -  6oo,0)  +  (kb  -  ftooK+1)^  p  ,n+i 


|n+l 


(6oo,0  +  ^ooP”+1) 


(A 


2(r0  +  KTpn+1) 


|n+l  p  |n+l 

rr 


ar+1) . 


(66) 


Substituting  the  update  equations  for  Px+1  and  Cx+1  in  (66)  we  obtain  an  expression 
for  Ex+1  in  terms  of  Ex,  Dx+1 ,  Cx  and  Px  given  by 


E(pn+1)  =  1  + 


At 


2  (t0  +  Krpn+1) 

+  (600, 0  +  KooPn+1)  + 


-1 


At((es,0  -  6oo,0)  +  X  ~  Kqo )pn+1) 
2  (r0  +  fivpn+1) 


2e, 


0 


n+l  _  -rY„n+l\-l 


Pai|n+1  =  P(F 


X<M  1- 


-1 


At 


2(r0  +  nTpn 


D  \n+l  -Cnx-  —Ex \n  -(l  +  - - At  Tp 

1  *  2e0  V  2(r0  +  fi)T^+1; 

P  in  +  |  At((ea,o  -  600, 0)  +  (Ks  -  /too)pra)  J  E  |n 


2(r0  +  KrP" 


(67) 
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where  Dx\n+1  is  calculated  from  (60).  Thus,  all  the  terms  on  the  right  side  are  known 
quantities,  and  we  can  use  equation  (67)  to  update  Ex.  Similarly  the  update  equation 
for  Ez  is 


E. 


|n+1  =  .F(pn+1)_1 
At 


D. 


in+i  -  cnz  -  °^ez r  -  [  i  + 


2e. 


o 


At 


2  (r0  +  Krpn+1 ) 


-l 


X 


1  - 


2(r0  +  nTpn 


Pz\n  + 


^^((^s, 0  ^oo,o)  T  {ks  Koo)P  ) 


E, 


(68) 


2(r0  +  Krpn) 

with  E{pn+l)  as  defined  in  (67).  As  before,  we  calculate  Dz\n+1  from  (62).  Hence  all 
terms  on  the  right  hand  side  are  known. 


4  The  forward  problem  for  a  Debye  medium:  First 
test  case 

As  discussed  in  Section  3.3  the  pressure  dependent  parameters  e*,  e^,  and  r*  are 
represented  as  a  mean  value  plus  a  perturbation  that  is  proportional  to  the  pressure  in 
equations  (47)  (as  discussed  in  [ABR02]).  In  our  first  test  case,  we  present  numerical 
simulations  for  a  Debye  medium  that  is  similar  to  water.  For  signal  bandwidths  in 
the  microwave  regime  the  dispersive  properties  of  pure  water  are  usually  modeled  by 
a  Debye  equation  having  a  single  molecular  relaxation  term.  The  mean  values  e*0, 
0,  and  Tq  for  this  Debye  model  [APM89]  are  given  by 

e*0  =80.1  (relative  static  permittivity), 

eoo  o  =5.5  (relative  high  frequency  permittivity), 

Tq  =  8.1  x  10”12  seconds, 

a*  =  1  x  10~5  mhos/meter, 

where  the  *  superscript  will  denote  the  true  values  of  all  the  corresponding  dielectric 
parameters.  Since  we  do  not  yet  have  experimental  data  to  determine  the  actual 
values,  we  instead  choose  trial  values  for  the  coefficients  of  pressure,  k*,  and  k*. 
As  a  first  approximation  we  take  these  trial  values  to  be  a  fraction  of  the  mean  values 
of  the  dielectric  parameters,  e*0,  e^Q,  and  Tq,  respectively.  Thus,  in  the  pressure 
region  we  choose 

=  0.6e:;0  =  48.06, 

<C  =  O.Se^o  =  4.4,  (70) 

k*  =  0.05tq  =  4.05  x  10-13. 

We  hope  to  determine  the  values  of  these  pressure  coefficients  from  appropriately 
designed  experiments  in  the  near  future.  We  have  assumed  that  the  effect  of  the 
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electric  field  on  the  acoustic  pressure  is  negligible.  We  are  interested  in  how  and 
to  what  extent  the  acoustic  pressure  can  change  the  reflected  electric  wave  from 
the  interface  at  z  =  z-i.  Moreover,  the  relaxation  parameter  is  a  characteristic  of 
the  material  in  [xq ,  z2\ ,  which  affects  the  transmission  of  electric  waves  through  this 
region.  Consequently,  we  will  examine  in  more  detail  how  To  can  effect  the  electric 
held.  The  electromagnetic  input  source  has  the  form 


J s(t,  x ,  z)  =  I(X1,X2)KZ)  sin (uc(t  -  3 10))  exp  (  - 

1 


t  —  3to 
to 


(71) 


to  — 


2vr  x  109 


;  uc  =  6tt  x  109  rad/sec;  fc  =  3.0  x  10J  Hz. 


The  Fourier  spectrum  of  this  pulse  has  even  symmetry  about  3.0  GHz.  The  com¬ 
putational  domain  is  defined  as  follows.  We  take  X0  =  (0,0.1),  Za  =  (0,0.15)  and 
Zb  =  (0.15,0.2).  The  number  of  nodes  along  the  2-axis  is  taken  to  be  320  and  the 
number  of  nodes  along  the  re-axis  is  taken  to  be  160.  The  spatial  step  size  in  both 
the  x  and  2  directions  is  Ax  =  Ax;  —  h  —  0.1/160.  From  the  CFL  condition  (55) 
with  the  Courant  number  r/cN  =  1/2  we  obtain  the  time  increment  to  be  At  zz  1.0417 
pico  seconds.  The  central  frequency  of  the  input  source  as  described  in  (71)  is  3.0 
GHz  and  based  on  the  speed  of  light  in  air,  c0  =  3  x  108  m/s,  we  calculate  the  corre¬ 
sponding  central  wavelength  to  be  Ac  =  (2nco)/u  =  0.1  meters.  The  antenna  is  half  a 
wavelength  long  and  is  placed  at  (aq,^)  x  zc,  with  2C  =  0,  aq  —  0.025  and  aq  =  0.075. 
We  use  PML  layers  that  are  half  a  wavelength  thick  on  all  fours  sides  of  the  compu¬ 
tational  domain  as  shown  in  Figure  3.  The  reflections  of  the  electromagnetic  pulse 
at  the  air-Debye  interface  and  from  the  acoustic  pressure  wave  are  recorded  at  the 
center  of  the  antenna  (xc,zc),  with  xc  =  0.05,  at  every  time  step.  This  data  will  be 
used  as  observations  for  the  parameter  identification  problem  to  be  presented  later. 
The  component  of  the  electric  held  that  is  of  interest  here  is  the  Ex  component.  Thus 
our  data  is  the  set 


E(q')  =  { EJvAt,  :rr. q'  j 


(72) 


The  windowed  acoustic  pressure  wave  as  defined  in  (49)  has  the  parameter  values 
Up  =  6.07T  x  105  rad/sec,  cp  =  1500  m/s,  and  thus,  \p  =  0.005.  The  location  of  the 
pressure  region  is  in  the  interval  (22,22  +  Ap)  =  (0.175,0.18). 

In  Figure  5  we  plot  the  electromagnetic  source  that  is  used  in  the  simulations  for 
the  Debye  model  (69).  We  plot  the  power  spectral  density  of  the  source  (71)  in  Figure 
6.  The  power  spectral  density  \Y\2  of  the  vector  Js  is  defined  to  be 


2  =  FFT(JS), 

|y|2  =  z  ■  z, 


(73) 
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Figure  5:  Forward  Simulation  for  a  Debye  medium  with  parameters  given  in  (69)  and 
(70).  The  source  as  defined  in  (71) 


where  FFT(JS)  is  the  fast  fourier  transform  of  the  vector  Js.  As  seen  in  Figure  6  the 
power  spectral  density  is  symmetric  about  3.0  GHz. 

In  Figures  7  and  8  we  plot  the  Ex  field  magnitude  at  the  center  of  the  antenna 
versus  time,  which  shows  the  electromagnetic  source,  the  reflection  off  the  air-Debye 
interface  and  the  reflection  from  the  region  containing  the  acoustic  pressure  wave. 
As  seen  in  these  plots  the  amplitude  of  the  reflection  from  the  acoustic  pressure  is 
many  orders  of  magnitude  smaller  than  that  of  the  initial  electromagnetic  source  as 
well  as  that  of  the  reflection  from  the  air-Debye  interface.  Thus,  it  is  an  interesting 
question  as  to  whether  the  reflections  from  the  acoustic  pressure  region  can  be  used 
for  identification  of  parameters. 

In  Figure  9  we  plot  the  power  spectral  density  of  the  data  in  Figure  7.  In  Fig¬ 
ures  10  -14  we  plot  the  Ex  field  magnitude  in  the  plane  containing  the  center  of  the 
antenna  versus  depth  along  the  z  axis.  Figures  10  and  11  shows  the  electromagnetic 
wave  penetrating  the  Debye  medium.  In  12  we  see  the  reflection  of  the  electromag¬ 
netic  wave  from  the  Debye  medium  moving  towards  the  antenna  and  the  Brillouin 
precursor  propagating  in  the  Debye  medium.  In  Figure  13  we  observe  the  reflection 
from  the  region  containing  the  acoustic  pressure  and  the  transmitted  part  of  the 
electromagnetic  source  travelling  into  the  absorbing  layer.  The  reflection  from  the 
acoustic  pressure  crosses  the  Debye  medium  in  Figure  14. 
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x  io3  The  power  spectral  density  of  the  Source 


Figure  6:  Fourier  transform  of  the  source.  The  transform  is  centered  around  the 
central  frequency  of  ui  =  6n  x  109. 


Figure  7:  Forward  simulation  for  a  Debye  medium  with  parameters  given  in  (69)  and 
(70).  Time  plot  of  the  Ex  component  of  the  electric  held.  This  is  the  data  that  is 
received  at  the  center  of  the  antenna.  This  plot  consists  of  the  original  signal  (source), 
the  reflection  off  the  Air-Debye  interface  and  the  reflection  due  to  the  pressure  wave. 
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Figure  8:  Forward  simulation  for  a  Debye  medium  with  parameters  given  in  (69) 
and  (70).  Time  plot  of  the  Ex  component  of  the  electric  field.  This  is  the  data  that 
is  received  at  the  center  of  the  antenna.  This  plot  consists  of  the  reflection  off  the 
Air-Debye  interface  and  the  reflection  due  to  the  pressure  wave. 


iq3  The  power  spectral  density  of  the  observed  Ex  field 


Figure  9:  Fourier  transform  of  the  data  in  Figure  7.  The  transform  is  centered  around 
the  central  frequency  of  ca  =  67T  x  109. 
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t=  7.29  xio-10  sec 


Figure  10:  Forward  simulation  for  a  Debye  medium  with  parameters  given  in  (69) 
and  (70).  Plot  of  the  magnitude  of  the  Ex  component  of  the  electric  field  against  z 
(meters).  The  input  source  is  propagating  towards  the  Debye  medium. 


t  =  1 .04  x  1 0'9  sec 


Figure  11:  Forward  simulation  for  a  Debye  medium  with  parameters  given  in  (69) 
and  (70).  Plot  of  the  magnitude  of  the  Ex  component  of  the  electric  held  against  z 
(meters).  The  input  signal  gets  partially  reflected  off  the  Debye  medium  and  partially 
transmitted  into  the  medium. 


t=  1 .25x  1 0“9  sec 


t=  1.46  x  1CT9  sec 


Figure  12:  Forward  simulation  for  a  Debye  medium  with  parameters  given  in  (69) 
and  (70).  Plots  of  the  magnitude  of  the  Ex  component  of  the  electric  field  against  z 
(meters).  In  both  the  plots,  the  reflection  off  the  air- Debye  interface  is  seen  moving 
towards  the  antenna,  while  the  transmitted  part  of  the  source  is  seen  propagating  in 
the  Debye  medium  towards  the  region  containing  the  pressure  wave. 
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Figure  13:  Forward  simulation  for  a  Debye  medium  with  parameters  given  in  (69) 
and  (70).  Plot  of  the  magnitude  of  the  Ex  component  of  the  electric  field  against  z 
(meters).  The  interaction  of  the  electromagnetic  wave  with  the  pressure  wave  causes 
the  electromagnetic  wave  to  partially  reflect  and  partially  transmit. 


4.1  Sensitivity  analysis 

Here  we  examine  the  system  dynamics  as  the  parameters  vary.  We  are  interested  in 
the  changes  produced  by  the  coefficents  of  pressure  in  the  polarization,  namely  the 
parameters  ks,  Kqo  and  kt.  In  Figure  15  we  plot  the  Ex  component  of  the  electric  field 
for  the  simulation  with  parameter  values  given  in  (69)  and  four  other  simulations  for 
which  all  the  parameters,  except  ks,  are  fixed  at  the  values  given  in  (69).  We  change 
the  value  of  ks  to  0,  0.3es,o,  0.6eSjo  and  0.9es,o,  repectively,  and  plot  the  electric  field 
values  for  each  of  the  corresponding  forward  simulations.  As  expected,  by  changing 
the  value  of  ns,  we  notice  a  change  in  the  magnitude  of  the  acoustic  reflection  observed 
at  the  center  of  the  antenna. 

In  order  to  determine  if  the  forward  problem  is  sensitive  to  changes  in  the  value  of 
the  pressure  coefficient  Kqo,  in  Figure  16  we  plot  the  absolute  value  of  the  differences 
in  the  Ex  field  component  between  simulations  in  which  all  the  parameters  except 
/too  are  fixed  at  the  values  given  in  (69).  The  solid  line  in  this  figure  represents 
the  difference  in  the  Ex  held  magnitude  between  the  forward  simulation  in  which 
Hoc  =  0.8600,0  and  the  simulation  in  which  Koo  =  0.  The  dashed  line  represents  the 
difference  in  the  Ex  held  magnitude  between  the  simulation  in  which  Koo  =  0.8eoo,o 
and  the  simulation  in  which  Koo  =  0.2600,0-  We  now  determine  if  the  forward  problem 
is  sensitive  to  changes  in  the  value  of  the  pressure  coefficient  nT.  In  Figure  17  we  plot 
the  absolute  value  of  the  differences  in  the  Ex  held  component  between  simulations  in 
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Figure  14:  Forward  simulation  for  a  Debye  medium  with  parameters  given  in  (69) 
and  (70).  Plot  of  the  magnitude  of  the  Ex  component  of  the  electric  field  against  z 
(meters).  The  reflection  from  the  acoustic  pressure  wave  impinges  on  the  air-Debye 
interface,  while  the  wave  transmitted  into  the  pressure  region  moves  into  the  right 
absorbing  layer. 


Figure  15:  Plot  of  the  magnitude  of  the  Ex  component  of  the  electric  field  against  z 
(meters)  for  different  values  of  the  parameter  ks.  The  other  parameters  are  fixed  at 
the  values  given  in  (69)  and  (70) 
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Figure  16:  Plot  of  the  absolute  value  of  the  difference  in  magnitude  of  the  Ex  compo¬ 
nent  of  the  electric  held,  between  simulations  in  which  Koo  =  0.8600,0  and  simulations 
with  different  values  of  against  z  (meters).  The  other  parameters  are  fixed  at  the 
values  given  in  (69)  and  (70). 


Timet  (secs)  xIO9 


Figure  17:  Plot  of  the  absolute  value  of  the  difference  in  magnitude  of  the  Ex  compo¬ 
nent  of  the  electric  held,  between  simulations  in  which  the  value  of  kt  =  0.05to  and 
simulations  with  different  values  of  kt,  against  z  (meters).  The  other  parameters  are 
fixed  at  their  true  values. 
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which  all  the  parameters  except  nT  are  fixed  at  the  values  given  in  (69).  The  solid  line 
in  this  figure  represents  the  difference  in  the  Ex  field  magnitude  between  the  forward 
simulation  in  which  nT  =  0.05to  and  the  simulation  in  which  kt  =  0.  The  dashed  line 
represents  the  difference  in  the  Ex  field  magnitude  between  the  simulation  in  which 
nT  =  0.05to  and  the  simulation  in  which  nT  =  0.2ro  From  these  plots  we  observe  that 
the  parameter  ks  seems  to  be  the  most  influential  in  the  wave  interaction,  whereas  the 
pressure  coefficients  Koo  and  nT  do  not  seem  to  be  as  influential.  This  observation  can 
be  supported  from  an  analysis  of  (40).  In  our  simulations  the  outgoing  and  reflected 
radiation  will  be  dominated  by  frequencies  near  the  center  frequency  3.0  GHz.  Thus, 
e*  will  be  dominated  by  frequencies  near  3.0  GHz.  In  this  problem  loTq  ~  0(  10  2) , 
and  e0  =  8  x  1CT12.  Rewriting  (40)  as 
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we  consider  each  of  the  three  terms  in  (74)  separately  to  determine  their  magnitude. 
The  magnitudes  of  the  corresponding  true  values  of  these  three  terms  (neglecting  the 
acoustic  pressure  terms)  are  approximately  given  as 
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Thus  we  see  that  e*  will  be  most  sensitive  to  the  static  permittivity  eSj0  and  the 
effects  of  Coop  and  r0  will  not  be  as  pronounced.  Also,  this  implies  that  e*  will  be 
more  sensitive  to  the  pressure  coefficient  ks  than  to  the  coefficients  Kqo  and  kt,  which 
coincides  with  the  observation  that  was  made  from  Figures  15,  16  and  17 

We  would  also  like  to  see  what  effect  the  acoustic  speed  and  frequency  have  on 
the  amplitude  of  the  reflection  from  the  acoustic  pressure  region.  In  Figure  18  we 
plot  the  acoustic  reflection  observed  at  the  center  of  the  antenna  for  different  values 
of  the  acoustic  frequency  and  in  Figure  19  we  plot  the  acoustic  reflection  for  different 
values  of  the  acoustic  speed.  By  changing  the  speed  or  the  frequency,  the  wavelength 
A p  changes  and  thus  the  size  of  the  interval  ( Z2 ,  £2  +  Ap)  in  which  the  pressure  wave  is 
generated.  The  amplitude  of  the  acoustic  reflection  appears  to  decrease  as  the  acoustic 
frequency  is  increased  and  increases  at  the  acoustic  speed  increases.  We  note  here 
that  our  windowed  pressure  wave  contains  only  one  wavelength  of  the  sinusoid.  Thus 
in  general  it  will  be  difficult  to  predict  precisely  how  the  reflections  will  behave  as  the 
speed  or  the  frequency  is  changed. 
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Figure  18:  Plot  of  the  reflection  from  the  pressure  region  that  is  measured  at  the 
center  of  the  antenna,  versus  time  t  (bottom),  corresponding  to  different  windowed 
acoustic  pressure  waves  (top)  of  varying  acoustic  frequencies.  The  other  parameters 
are  fixed  at  the  values  given  in  (69)  and  (70).  We  note  that  the  amplitude  of  the 
reflection  changes  as  the  frequency  of  the  pressure  wave  changes. 
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Figure  19:  Plot  of  the  reflection  from  the  acoustic  region  that  is  measured  at  the 
center  of  the  antenna,  against  time  t  (bottom),  corresponding  to  different  windowed 
acoustic  pressure  waves  (top)  of  varying  acoustic  speeds.  The  other  parameters  are 
fixed  at  the  values  given  in  (69)  and  (70).  We  observe  that  the  magnitude  of  the 
reflection  is  affected  by  the  speed  of  the  acoustic  wave. 
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5  Parameter  estimation  and  statistical  inference: 
The  inverse  problem 

In  our  forward  simulations  the  signal  that  we  record  at  the  center  of  the  antenna  is  a 
set  of  measurements  of  Ex,  the  x  component  of  the  electric  held,  at  the  point  (0 ,zc) 
in  our  computational  domain  at  discrete,  uniformly  spaced  intervals  of  time.  This 
signal  has  components  consisting  of  the  input  source,  the  reflection  at  the  air-Debye 
interface  and  another  reflection  from  the  region  that  contains  the  pressure  wave  as 
seen  in  Section  4.  The  signal  is  a  function  of  the  various  dielectric  properties  of  the 
Debye  medium;  namely,  the  static  relative  permittivity  eS}o,  the  infinite  frequency 
permittivity  e^o,  the  relaxation  time  To,  the  conductivity  a,  and  the  three  pressure 
coefficients  ks,  and  kt.  We  collect  all  these  quantities  together  to  define  the 
parameter  vector 

Q  [^s,0>  £oc,0;  ^si  ^ooi  ^t\  ■  (76) 

Thus,  the  signal  recorded  in  the  forward  simulation  is  a  function  of  q,  and  changing 
the  value  of  q  will  change  the  signal  that  is  recorded.  In  general,  such  signals  are 
usually  obtained  in  an  experimental  setting  conducted  in  a  laboratory  with  physical 
equipment  such  as  electric  pulse  generators  to  generate  electric  pulses,  piezoelectric 
transducers  that  generate  the  acoustic  pressure  waves  and  antennas/receivers  that 
record  the  electric  held  intensities.  Moreover,  such  signals  generated  in  an  experi¬ 
mental  setting  are  usually  noisy,  with  noise  arising  through  the  equipment  that  is 
used;  namely  the  resistor,  antenna/receiver  and  the  transducer,  in  our  case.  Implicit 
in  such  collection  of  measurements  is  the  assumption  that  there  exist  true  values  of 
all  the  parameters  that  characterize  the  medium  to  be  interrogated.  We  will  denote 
the  corresponding  true  parameter  vector  by  q*. 

In  order  to  simulate  an  experimental  situation  (i.e.,  generate  typical  “data”)  we 
assume  that  there  exist  true  values  of  all  the  parameters  involved  in  our  problem, 
and  we  use  these  true  values  (which  of  course  would  not  be  known  in  an  experimental 
setting!)  in  our  forward  simulation,  viz.  the  FDTD  method  and  record  the  signal 
observed  at  the  antenna  as  described  above.  We  then  add  noise  to  the  generated  signal 
to  produce  a  noisy  signal  which  will  form  the  observations  or  data  that  we  will  use  in  an 
inverse  problem,  as  a  substitute  for  data  generated  in  the  laboratory.  The  goal  of  our 
electromagnetic  interrogation  technique  is  to  identify  or  estimate  the  true  parameter 
q*  that  characterizes  the  particular  Debye  medium  being  interrogated,  from  data 
that  consists  of  a  signal  recorded  at  the  center  of  the  antenna.  The  reason  for  using 
simulated  data  in  the  inverse  problem  is  to  validate  the  methods  first  on  “data”  from 
known  parameters  in  a  setting  with  known  noise.  If  we  cannot  estimate  the  dielectric 
parameters  from  data  that  is  constructed  via  numerical  simulations  of  our  discrete 
model,  then  the  reconstruction  of  these  parameters  from  actual  experimental  data 
usually  will  not  be  feasible. 

We  can  state  the  inverse  or  parameter  estimation  problem  that  we  will  attempt  to 
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solve  as  follows:  Using  observations  or  data  (electric  field  intensities  containing  noise 
collected  at  the  center  of  the  antenna),  determine  an  estimate  q,  of  the  true  parameter 
q*,  that  belongs  to  an  admissible  set  Q,  so  that  the  solution  of  the  forward  problem 
with  q  =  q  best  describes  the  data  that  is  collected.  When  we  have  solved  this 
deterministic  problem  (as  will  be  detailed  below),  it  will  also  be  necessary  to  specify 
reliablility  of  our  estimates,  i.e.,  can  we  specify  measures  of  uncertainty  related  to 
the  estimate  q  of  q*?  It  is  important  to  note  that  the  uncertainty  in  consideration  is 
inherent  in  the  method  of  producing  the  estimates  as  well  as  in  the  process  of  data 
collection.  Such  measures  of  uncertainty  will  be  specified  by  means  of  confidence 
intervals.  These  intervals  are  a  probability  statement  about  the  procedure  which 
is  used  to  construct  estimates  of  the  parameters.  Thus  we  are  led  in  a  completely 
natural  way  to  stochastic  or  probabilistic  aspects  of  estimates  from  a  deterministic 
problem  solved  with  deterministic  algorithms  [BB99] .  The  statistical  error  analysis 
that  we  will  present  here  is  based  on  standard  statistical  formulations  as  given  in 
[DG95], 

We  first  consider  the  two  different  ways  in  which  we  can  add  noise  to  our  signal, 
i.e.,  the  values  of  the  x  component,  Ex,  of  the  electric  field  observed  at  the  center  of 
the  antenna,  (0,  zc).  for  different  times  tk  =  kAt ,  k  —  0, 1, ... ,  M .  We  represent  this 
signal  as  the  vector 


(i)  E(q*)  =  {E*k}^=1  =  (Ex(ti,  0,  zc;  q*), . . . ,  Ex{tM,  0,  ^c;  q*))'  , 


n 


*\T 


(77) 


Q  ^oo,0 ?  ^~0  ?  ®  ^t)  5 


q*  being  the  true  parameter  vector. 


1.  Relative  random  noise  (RN)  :  In  this  case  the  amplitude  of  the  noise  that  is 
added  is  proportional  to  the  size  of  the  signal,  {Efi)^f=v  Our  simulated  data  is 

Ol=E-t(l  +  v,fk),  k=l,...M,  (78) 

where  fj  =  {r]k}^=\  are  independent  normally  distributed  random  variables  with 
mean  zero  and  variance  one,  i.e.,  ~  J\f( 0,  l),k  =  1, . . .  ,M.  We  express  the 
relative  magnitude  of  the  noise  as  a  percentage  of  the  magnitude  of  E(q*)  by 
taking  two  times  the  value  of  u  as  the  size  of  the  random  variable.  For  example, 
v  =  0.005  corresponds  to  1%  relative  noise  [BBL00].  This  noise  model  does  not 
produce  constant  variance  across  the  samples. 

2.  Constant  variance  noise-.  Since  constant  variance  is  most  conveniently  assumed 
in  standard  error  analysis,  we  will  consider  estimates  obtained  from  an  inverse 
problem  applied  to  simulated  data  which  contains  constant  variance  random 
noise.  The  data  that  we  try  to  fit  in  this  case  is 


Ot  =  E*k  +  k  =  1, . . .  M, 


(79) 
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where  as  in  (78),  rjk  ~  A/”(0, 1),  k  =  1, . . . ,  M.  The  constant  /3  is  taken  to  be 
the  product  va,  where  a  is  a  scaling  factor  which  is  chosen  so  that  the  signal 
to  noise  ratio  (SNR)  of  the  data  defined  in  (79)  corresponds  to  data  defined  in 
(78)  with  noise  level  v.  This  justifies  comparing  results  obtained  with  the  two 
different  ways  of  adding  comparable  noise  to  our  signal.  This  will  be  addressed 
in  more  detail  in  Section  5. 

The  vector  Os  =  will  be  our  data  or  observations,  and  the  noise  fj  =  {rjk}kLi 

will  be  referred  to  as  the  measurement  errors  corresponding  to  the  observations. 
The  statistical  error  analysis  that  we  now  develop  is  only  applicable  to  the  case  of 
constant  variance  noise.  Thus,  for  the  rest  of  this  section,  we  will  assume  that  our 
observations  are  of  the  form  (79).  The  sample  observations  Os,  which  are  in  general 
obtained  from  experiment,  are  a  realization  of  the  corresponding  random  variable 
O  =  (Oi,02,...,Om)t  which  can  be  seen  to  be  a  transformation  of  the  random 
variable  fj.  Thus 

Ok  =  E*k  +  / 3r]k ,  k  =  1, . . .  M,  (80) 

is  a  stochastic  process  and  has  a  multivariate  normal  distribution  with  mean  vector 
E(q*)  defined  in  (77,  i)  and  covariance  matrix  (32Im-  The  statistical  model 

O  ~  AV{E(q*), /32Zm},  (81) 

is  a  formal  representation  of  the  population  of  all  possible  realizations  of  0\,  02,  — ,  Om 
that  can  be  observed.  When  we  collect  data,  we  are  observing  a  single  realization 
of  Of, . . . ,  OsM  i.e.,  a  sample.  Our  objective  then  is  to  estimate  the  true  value  of  the 
parameter  q*  by  collecting  data,  i.e.,  observing  a  single  realization  of  O  as  well  as 
accounting  for  the  fact  that  a  different  realization  will  produce  a  different  estimate  of 
q*.  We  would  like  to  learn  about  the  true  value  q*  (which  determines  the  nature  of 
the  population)  from  a  sample,  as  well  as  indicate  the  certainty  (or  uncertainty)  that 
we  can  associate  to  our  knowledge  of  this  parameter.  This  process  of  making  state¬ 
ments  about  a  population  of  interest  on  the  basis  of  a  sample  from  the  population  is 
called  statistical  inference. 

The  unknown  true  parameters  vector  q*  will  be  estimated  by  means  of  a  suit¬ 
able  function  q(O)  of  the  observations  O.  The  function  q(O)  is  a  random  vari¬ 
able  and  is  referred  to  as  an  estimator.  When  evaluated  at  a  particular  realization, 
01,01,...  ,OsM,  this  estimator  yields  a  numerical  value  that  gives  information  about 
the  true  value  of  the  parameter  q*.  For  a  particular  realization  Os  of  the  random 
variable  of  observations  O,  we  call  q(Os)  an  estimate.  In  this  report  we  will  consider 
the  method  of  least  squares  to  obtain  estimates  for  the  parameters  and  to  calculate 
confidence  intervals  for  the  estimates  by  linearizing  around  the  estimate  q(Os).  In 
this  case  q(O)  will  be  the  least  squares  estimator  q0Ls(0).  Associated  with  the  ran¬ 
dom  variable  qoLs(O)  is  the  probability  space  Q  =  ( Q ,  B ,  m),  where  Q  is  the  sample 
set  of  all  admissible  parameters  q,  B  is  the  u-algebra  of  events  and  m  is  the  probabil¬ 
ity  measure  (or  distribution).  Thus  we  look  for  the  least  squares  estimator  qoLs(O) 
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on  0  such  that 


qoLs(O)  =  arg  min  J(q), 

q  on  Q 

l  M 

J(q)  =  2  ^  ^  I Ex(tk,  0,  zc\  q)  Ok 

z  k= 1 


(82) 


where  E(q)  =  (Ex(ti,  0,  zc;  q),  Ex(t2,  0,  zc;  q), . . . ,  Ex(tM,  0,  zc;  q))T  will  be  generated 
by  our  forward  simulation.  We  will  refer  to  the  vector  E(q)  as  our  simulations.  For  a 
particular  realization  Os  of  the  random  variable  O  of  observations,  this  process  will 
yield  the  least  squares  estimate  qoLs(Os)  G  Q,  where 


qoLs(Os)  =  arg  min  J5(q), 
qe  Q 

1  M 

j5(q)  =  2  S  \ESk,  o,  ^c;  q)  -  o*k |2 , 

Z  k= 1 


(83) 


5.1  Calculation  of  confidence  intervals  via  linearization 

In  general,  the  mappings  between  the  parameters  q  and  the  simulations  E  are  non¬ 
linear.  Let  Ek  =  Ex(tk,  0,  zc\  q).  The  functions 

Ei{qi, . . .  ,  qi)  —  Ex(ti,  0,  zc;  q), 

E2(qu  •  •  •  ,Qi)  =  Ex(t2,  0,  zc;  q), 


EM(qi,  ...,qi)  =  Ex(tM,  0,  zc;  q), 


denote  real- valued  differentiable  functions  of  the  unknown  parameters  q  =  (qu  . . . ,  qi)T , 
where  1  <  l  <  7  depending  on  how  many  parameters  we  attempt  to  identify.  Let 
q  =  q*  +  Dq,  where  the  corrector  term  Z>q  =  (Agi, . . .  A qi)T  is  unknown.  Linearizing 
the  functions  in  (84)  around  the  true  parameter  q*  by  using  the  Taylor  expansion  we 
obtain 


Ek(qi,  ...qi)  =  Ek(ql  +  Aqu...,  q*  +  A  ®) 

1 


=  Ek(  q* 


dEk 

hi  dqi 


Am,  k  —  1 . . .  M. 


Let  us  define 


and 


Vc  =  (Oi  -  Ex{ q*), .  .  .  ,0M  -  Em( q*))T, 


T(q*) 


/ 


8Ei 

dqi  q* 


l  9Em 

\  dqi  q* 


(85) 


(86) 


(87) 
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Thus,  in  order  to  estimate  the  unknown  corrector  terms  by  the  method  of  least 
squares,  the  estimator  Pols  is  determined  by 

Pols  =  min \{Ve  -  A’(q*)Pq)T(Pe  -  *(q*)Pq),  (88) 

which  is  the  linearized  version  of  (82).  This  involves  solving 


9  (PeTPe  -  2PjT(q*)P0LS  +  V^LSX(q*)TX(q*)T>OLS) 

— ^ - - - J-  =  0, 

<9Pols 

and  the  corresponding  optimality  condition  yields  the  estimator 


(89) 


Pols  =  (T’(q*)TA’(q*))-1T(q*)TPc.  (90) 

We  note  that  the  expected  value  of  Pols  is  Exp(Pols)  =  0,  and  the  covariance  matrix 
of  Pols  and  hence  also  of  qoLS  is 

Cov(Pols)  =  (T’(q*)TT’(q*))-1T(q*)TCov(Pe)T’(q*)(T(q*)TT’(q*))-1 

= /?2(T’(q*)TT(q*))'1  =  Cov(qoLs), 


where  qoLS  =  q*  +  Pols,  and  (3  is  defined  in  (79).  Thus  the  least  squares  estimator 
is  a  random  variable  with  a  multivariate  normal  distribution  which  after  linearization 
can  be  approximately  represented  as 

qoLs(O)  ~  Mi  (q*,  /32 [T’T(q*)T’(q*)]~1)  .  (92) 

However,  when  data  is  generated  in  an  experimental  setting,  we  do  not  have  any 
knowledge  of  the  true  parameter  vector  q*  and  hence  we  cannot  calculate  (X (q*)TT (q*)) 
Thus,  we  further  approximate  our  least  squares  estimator  as  having  the  multivariate 
normal  distribution 

qOLs(O)  ~  Ml  ^q0Ls(Of>),/3oLS[^7(qOLs(0S))^(qOLs(C)6))]^1)  ,  (93) 

which  is  obtained  by  repeating  the  linearization  process  around  the  estimate  qoLs(Os) 
obtained  from  a  realization  Os  of  the  random  variable  O  of  observations.  The  value 
of  /3qls  ch°sen  as 


Pols  ~ 


^  I Ex(tk,  0,  2;c;  qoLs(Os))  —  Osk 


which  is  the  minimum  value  of  the  least  squares  objective  function  scaled  by  2/(M— /). 

Whenever  an  estimate  based  on  data  is  reported,  it  should  be  accompanied  by 
an  assessment  of  uncertainty  based  on  the  sampling  distribution.  To  this  end  we 
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will  construct  confidence  intervals  for  all  the  estimates  that  we  will  provide.  The 
approach  we  use  is  to  look  at  the  standard  error  (SE)  for  each  of  the  l  components 
of  the  estimate  qoLs(Os),  which  is  given  for  its  jth  component  by  the  jth  diagonal 
term  of  the  covariance  matrix,  i.e., 


SE(goLSj)  =  \/var(g0LSj(Os)) 

v _  (95) 

=  \J Pols  ((;t’(qoLs(Os))rT(q0Ls(Os)))-1)ii  j  =  l,...,l. 

We  construct  the  intervals 

CIj  =  (?ols,j(Os)  —  l-96S'£,(g0LSj),  Qohs,j(Os)  +  1.96SE(qOLS,j)) ,  j  —  1, . . . 

(96) 

for  which 

m  ({<Zolsj(Os)  —  1.96SE(q0Ls,j)  <  qoLS,j(0)  <  Qohs,j(Os)  +  1.96SE(qoLs,j)})  =  0.95, 
j  —  1, . . .  ,1, 

(97) 

where  m  is  the  probability  measure  for  the  probability  space  Q,  and  M  is  sufficiently 
large  that  one  can  use  a  Gaussian  A/”(0, 1)  distribution  in  computing  confidence  in¬ 
tervals.  These  intervals  CL,-  are  called  confidence  intervals  with  confidence  level  0.95 
or  the  95%  confidence  interval. 

Remark  2  The  confidence  intervals  are  a  probability  statement  about  the  procedure 
by  which  an  estimate  is  constructed  from  a  sample  of  the  population.  They  are  not  a 
probability  statement  about  the  true  parameter  q*  which  is  fixed.  If  we  could  construct 
the  confidence  intervals,  from  our  estimation  procedure,  for  all  possible  data  samples 
of  size  M,  then  95%  of  such  intervals  would  cover  the  true  parameter  values  q*. 
However,  for  a  particular  data  sample  Os,  the  confidence  intervals  constructed  as 
above  may  or  may  not  cover  the  true  parameter  q* .  What  we  can  state  is  that  we  are 
95%  confident  that  the  confidence  intervals  constructed  by  our  estimation  procedure 
would  cover  q* . 


6  Parameter  estimation:  First  test  problem 


We  now  attempt  to  estimate  the  dielectric  parameters  for  the  Debye  medium  defined 
in  (69)  and  (70).  We  restate  the  values  of  all  the  parameters  below 


f* 

fcs,0 

=  80.1 

K*s 

=  48.06 

f* 

too,0 

=  5.5 

n* 

=  4.4 

ro 

=  8.1  x  10”12 

OO 

<7* 

=  1  x  10“5 

< 

=  4.05  x  10 

(98) 
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We  will  refer  to  the  values  (98)  as  the  true  values  of  the  parameters.  We  will  solve 
the  corresponding  least  squares  optimization  problem  using  two  different  methods; 
the  Nelder-Mead  algorithm,  which  is  a  simplex  based,  gradient  free  method,  and 
the  Levenberg-Marquardt  trust  region  method  that  uses  forward  finite  differences  to 
calculate  the  gradient.  From  our  estimates  for  all  the  seven  parameters  we  will  show 
that  our  problem  is  sensitive  to  only  three  of  the  seven  parameters.  We  will  then 
try  to  estimate  two  or  three  of  the  parameters  to  which  our  problem  is  sensitive. 
The  particular  implementation  of  the  Levenberg-Marquardt  and  the  Nelder-Mead 
algorithm  that  are  used  in  this  report  are  based  on  the  corresponding  algorithms 
presented  in  [Kel99] 

6.1  Simulation  results:  The  Nelder-Mead  method 

We  first  present  parameter  estimation  results  using  the  Nelder-Mead  simplex  based 
algorithm  which  is  used  for  optimization  in  noisy  problems.  The  Nelder-Mead  simplex 
algorithm  keeps  a  simplex  S  of  approximations  to  an  optimal  point.  In  this  algorithm 
the  vertices  of  the  simplex  where  /  is  the  size  of  q,  are  sorted  according  to 

the  least  squares  objective  function  values 


Js(qi)  <  Js(q2)  <  •  •  •  Js(qm)' 


(99) 


q!  is  called  the  best  vertex  and  q;+1  the  worst.  The  algorithm  attempts  to  replace 
the  worst  vertex  q;+1  with  a  new  point  of  the  form 


q(^nm)  =  (1  +  ^nm  )  q  Vnm  qj+ 1 5 
where  q  is  the  centroid  of  the  convex  hull  of  {qj}[=i 


i= 1 


(100) 


(101) 


The  value  of  unm  is  selected  from  a  sequence  of  values  that  are  computed  by  performing 
several  different  operations  on  the  simplex.  In  general,  the  Nelder-Mead  algorithm  is 
not  guaranteed  to  converge,  even  for  smooth  problems.  The  failure  mode  is  stagnation 
at  a  nonoptimal  point.  For  further  details  we  refer  the  reader  to  [Kel99] . 

In  the  first  test  we  will  try  to  estimate  all  of  seven  parameters  in  q*.  The  initial 
simplex  is  chosen  to  have  values  that  are  5-10%  larger  than  the  true  parameter  values 
given  in  (98).  We  do  not  add  any  noise  to  our  data  that  is  used  in  the  parameter 
estimation  in  this  first  test.  The  final  estimates  from  the  Nelder-Mead  algorithm 
are  presented  in  Table  1.  The  details  of  the  simulation  are  plotted  in  Figures  21- 
22.  Figure  21  plots  the  various  features  of  the  Nelder-Mead  run  as  functions  of  the 
number  of  iterations.  The  algorithm  is  terminated  when  the  maximum  least  squares 
function  value  difference  between  any  two  points  in  the  simplex  is  less  than  10~8. 
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Parameter 

True  Values 

Final  Estimates 

0,0 

80.1 

80.1111 

^oo,0 

5.5 

5.7461 

P) 

8.1  x  10~12 

8.1283  x  10~12 

a 

1.0  x  10"5 

1.1972  x  10'5 

Ks 

48.06 

48.1227 

^oo 

4.4 

4.9148 

Kt 

4.05  x  10”13 

4.1935  x  10~13 

Table  1:  Final  estimates  for  all  seven  param¬ 
eters  for  water  using  Nelder-Mead.  The  ini¬ 
tial  parameter  simplex  has  values  that  are 
5-10%  larger  than  the  true  values. 


Function  Difference 

7.6979  x  10~9 

Max  Oriented  Length 

2.4986  x  10“3 

Average  least  squares  value 

9.8451  x  10“6 

Best  Least  squares  value 

9.8392  x  10“6 

L2  norm  of  the  gradient 

1.008  x  10"2 

L°°  norm  of  gradient 

8.756  x  10“3 

Number  of  forward  solves 

453 

Table  2:  Status  of  the  Nelder-Mead  Simula¬ 
tion  at  the  283rd  iteration. 
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O 
c n 


<to 


Ks  =  48.06  (eg)  x  10  12x0=  8.1e-12 


Figure  20:  Variation  of  eS)o,  To  and  ks  with  respect  to  iteration  number.  Out  of  all 
seven  parameters  these  three  parameters  are  relatively  well  estimated. 


The  algorithm  converges  in  283  iterations.  Table  2  presents  the  final  status  of  the 
simulation.  In  Figure  21,  the  maximum  oriented  length  is  defined  as 

Max  Oriented  length  =  max  llqi  —  qJI.  (102) 

2<i<m 

Also  the  gradient  in  this  case  refers  to  the  simplex  gradient.  As  can  be  seen  from 
Figure  22  the  parameters  e^o,  cr  and  the  two  pressure  coefficients  Kqo  of  e^o  and 
nT  of  To  are  difficult  for  the  inverse  problem  to  identify.  On  the  other  hand,  from 
Figure  20  we  see  that  the  parameters  eSio  and  its  pressure  coefficient  ks  as  well  as 
the  parameter  To  are  relatively  well  estimated.  Since  the  inverse  problem  is  only  able 
to  identify  three  parameters  in  the  absence  of  noise  we  do  not  expect  to  see  any 
improvement  when  noise  is  added  to  the  data!  Thus  we  will  henceforth  concentrate 
on  the  identification  of  es,o,  and  ks. 

We  will  first  attempt  the  identification  of  eSio  and  ks.  In  a  first  test,  we  fix  all 
the  other  5  parameters  (including  r0)  at  the  true  values  and  attempt  to  identify  eSi0 
and  ks.  Table  3  presents  the  results  for  this  test.  As  can  be  seen  from  Table  3,  with 
the  other  5  parameters  fixed  at  their  true  values  we  do  not  have  any  difficulty  in 
identifying  eS)o  and  ks  accurately.  In  Table  3,  FC  denotes  the  function  count,  i.e.,  the 
number  of  times  the  least  squares  objective  function  is  evaluated.  In  Figure  23  we 
plot  the  values  of  eSj0  (left)  and  ks  (right)  over  all  iterations,  for  the  three  different 
initial  simplicies  with  values  that  are  5%-10%  lower  than  q*  (o-),  5%  -10%  higher 
than  q*(-  -),  and  50%-60%  lower  than  q*(-). 

In  Figure  24  we  plot  the  details  of  the  Nelder-Mead  simulations  for  the  three 
different  initial  simplices  tabulated  in  Table  3.  In  each  case  the  algorithm  converges 
when  the  difference  between  function  values  is  less  than  10~8.  For  the  second  case 
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Function  Differences 


5 


Least  Squares  Function  Value 


.2 


L°°  Norm  of  Gradient 


.4 


Iterations 


Figure  21:  Nelder-Mead  for  the  estir 


Max  Oriented  Length 


Average  Least  Squares  value 


Iterations 


of  seven  parameters  for  water. 


k  =  4.4  k  =  4.05e-13 

xlCf13  1 


Iterations  Iterations 


Figure  22:  Variation  of  e^o,  <r  and  estimates  of  the  two  pressure  coefficients,  namely, 
/too  and  kT,  with  respect  to  iteration  number.  As  can  be  seen  our  inverse  problem  is 
not  sensitive  to  these  four  parameters  and  their  identification  is  difficult  even  in  the 
absence  of  noise. 
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e  =  80.1 

s,0 


k*  =  48.06 

s 


Figure  23:  Plots  of  eS)0  (left)  and  ks  (right)  against  the  number  of  iterations  for 
simulations  in  which  the  other  five  parameters  are  fixed  at  their  true  values  and  no 
noise  is  added  to  the  data,  (o-)  0.95q*,  (-  -)  1.05q*,  (-)  0.5q*. 


Table  3:  Parameter  estimation  of  es  0  and  ks  for  different  initial  simplices  with  no 
added  noise. 


qo 

Iter 

es,0 

ks 

|J^i  -  Jo  1 

Jy 

l|VP'||L2 

FC 

0.95q* 

33 

80.0999 

48.0607 

4.7019xl0-9 

1.4560x  10-8 

0.0204 

65 

1.05q* 

32 

80.1 

48.06 

9.5079xl0-9 

7.1205xl0-27 

0.0096 

66 

0.5q* 

59 

80.1001 

48.0603 

8.4746xl0-9 

5.5171xl0-9 

0.0130 

112 
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where  the  initial  simplex  has  points  with  5%-10%  larger  values  than  the  true  values, 
we  obtain  an  unusually  small  least  squares  function  value  of  the  order  of  10-2'.  We 
believe  this  to  be  an  accidental  artifact  of  the  algorithm.  In  Figure  25  and  26  we  plot 
the  least  squares  objective  function  for  values  of  es  0  and  ks  that  differ  by  0%-10%  from 
the  true  values.  The  other  five  parameters  are  kept  fixed  at  their  true  values.  From 
Figures  25  and  26,  we  determine  that  our  inverse  problem  is  more  sensitive  to  eS)o  than 
to  the  pressure  coefficient  ks.  This  is  expected  as  the  static  relative  permittivity,  eSio 
is  more  influential  in  the  system  dynamics  than  the  pressure  coefficients,  and  hence 
is  also  more  important  in  identifying  and  characterizing  the  material. 

As  explained  in  Section  5,  in  a  more  realistic  situation  we  will  not  have  access 
to  the  true  values  of  any  of  the  parameters.  Hence,  we  now  fix  five  of  the  seven 
parameters  at  values  that  have  a  relative  error  from  the  true  values  of  about  10%- 
50%.  We  use  the  values 

eoo,o  —  6.0, 

r0  =  10.0  x  10~12, 

a  =  1.5  x  10~5,  (103) 

/too  =  4.8, 

Kr  =  5.0  X  10-13, 

in  our  simulations.  Such  choices  for  the  parameters  are  typical  in  testing  algorithms 

[BK89].  We  use  different  initial  guesses  for  eS;o  and  ks.  Also  in  an  experimental  setting 
noise  is  introduced  into  the  data,  as  discussed  earlier.  To  simulate  the  experimental 
process  in  data  collection  we  first  add  relative  random  noise  to  our  data,  and  we  use 
the  fixed  values  (103)  in  the  inverse  problem.  The  results  are  tabulated  in  4.  From 
Table  4  we  see  that  using  the  inverse  problem  we  are  able  to  estimate  eSjo  but  we 
are  unable  to  reconstruct  ks  with  good  accuracy.  In  Figure  27  we  plot  the  values  of 
eSto  (left)  and  ks  (right)  over  all  iterations,  In  Figure  28  we  plot  the  details  of  the 
Nelder-Mead  simulation  for  the  results  in  Figure  27  and  Table  4  which  show  the  final 
estimates  of  eSio  and  ks  for  water.  The  other  five  parameters  are  fixed  at  the  values 
given  in  (103)  and  relative  noise  of  varying  levels  is  added  to  the  data. 

To  understand  why  ks  is  not  recovered,  we  note  that  the  fixed  value  of  r0  in  (103) 
has  a  relative  error  of  23%  from  its  true  value.  Since  the  reflections  from  the  acoustic 


Table  4:  Parameter  estimation  of  eSto  and  ks  for  0%,  1%,  3%  and  5% 
relative  noise. 


%  RN 

Iter 

£s,0 

re+i-jg] 

J6 

I|vj6’||l2 

FC 

0.0 

37 

80.3636 

66.0486 

8.5719xl0-9 

2.5594 

0.0076 

73 

1.0 

41 

80.3638 

66.0422 

3.2939xl0-9 

25.2195 

0.0098 

80 

3.0 

40 

80.3640 

66.0280 

8.6023xl0-9 

208.1146 

0.0080 

77 

5.0 

40 

80.3643 

66.0157 

6.3563xl0-9 

574.4429 

0.0060 

79 
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Iterations 


Iterations 


Figure  24:  Nelder-Mead  for  the  estimation  of  eS)o  arid  ks  for  water.  The  other  five 
parameters  are  fixed  at  their  true  values  and  no  noise  is  added  to  the  data,  (o-) 
0.95q*,  (-  -)  1.05q*,  (-)  0.5q*. 
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Least  Squares  Objective  Function 


Figure  25:  The  Least  Squares  Objective  Function  for  different  values  of  eSio  and  ks. 
We  can  see  from  this  plot  that  the  problem  is  more  sensitive  to  eS)o  than  to  the 
pressure  coefficient  ks. 


Least  Squares  Objective  Function  Values 


Figure  26:  Cross-sections  of  the  Least  Squares  Objective  function  versus  ks  (left)  and 
versus  eSjo  (right).  Note  the  difference  in  scales  between  the  two  figures. 
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e*  =  80.1  k*  =  48.06 

s,0  s 


Figure  27:  Plots  of  <%0  (left)  and  ks  (right)  against  the  number  of  iterations  for 
simulations  in  which  the  other  five  parameters  are  fixed  at  values  that  have  10%-50% 
relative  error  from  the  true  values  (103)  and  relative  noise  of  varying  levels  is  added 
to  the  data. 


region  depend  strongly  on  the  value  of  r0,  a  relative  error  of  23%  is  too  high  to  get 
a  good  estimate  of  the  pressure  coefficients.  To  test  this  we  lower  the  value  of  r0  to 
8.91  x  10-12.  So  instead  of  the  values  in  (103)  we  will  use  the  values 

^oo,o  6.0, 

r0  =  8.91  x  10~12, 

O-  =  1.5  x  10"5,  (104) 

Kqq  4.8, 

KT  =  5.0  X  10-13. 

This  value  of  r0  has  a  relative  error  of  10%  from  its  true  value  Tq.  With  this  new 

fixed  value  for  To  we  attempt  the  inverse  problem  again  to  reconstruct  eSio  and  ks. 
These  results  are  presented  in  Table  5.  From  Table  5  we  observe  that  the  estimated 
value  of  ks  is  closer  to  the  true  value  than  in  the  previous  test.  Thus  we  can  obtain 
accurate  estimates  of  ks  by  choosing  values  of  r0  that  are  closer  to  its  true  value.  We 
can  conclude  that  if  the  value  of  tq  is  not  accurately  known,  then  the  identification 
of  ks  is  difficult.  The  identification  of  eS;o  on  the  other  hand  does  not  appear  to  be 
sensitive  to  the  fixed  value  of  r0  that  is  chosen.  The  combined  identification  of  both 
e^o  and  ks  appears  to  be  quite  insensitive  to  the  level  of  relative  noise  that  is  added 
to  the  data.  The  results  of  the  Nelder-Mead  run  and  the  variation  of  the  parameters 
over  all  iterations  are  plotted  in  Figures  29  and  30,  respectively. 
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Function  Differences  Max  Oriented  Lenqth 


Iterations  Iterations 

Figure  28:  Nelder-Mead  for  the  estimation  of  eSio  and  ks  for  water.  The  other  five 
parameters  are  fixed  at  the  values  given  in  (103)  and  relative  noise  of  varying  levels 
is  added  to  the  data.  0%  (-),  1%  (-  -)  3%  (o-)  and  5%  (+-) 
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Table  5:  Parameter  estimation  of  eS)0  and  ks  for  0%  -  10%  added  relative 
random  noise.  The  values  of  the  other  five  parameters  are  fixed  at  values 
given  in  (104). 


%  RN 

Iter 

ks 

ij;ui  -  yi 

J6 

l|VJs||i» 

FC 

0.0 

38 

80.2165 

55.2945 

1.5487xl0-9 

0.4369 

0.0013 

73 

0.1 

38 

80.2165 

55.2949 

7.2106xl0-9 

0.6551 

0.006310 

75 

0.3 

39 

80.2166 

55.2946 

5.4008xl0-9 

2.4674 

0.005495 

76 

0.5 

39 

80.2167 

55.2939 

1.8372xl0-9 

6.1140 

0.009746 

76 

1.0 

39 

80.2169 

55.2914 

4.5901xl0-9 

23.2556 

0.0075 

76 

3.0 

41 

80.2173 

55.2826 

4.2341xl0-9 

206.4680 

0.0087 

78 

5.0 

40 

80.2178 

55.2729 

8.2103xl0-9 

573.1135 

0.0079 

78 

10.0 

39 

80.2191 

55.2509 

2.5361xl0-9 

2292.24730 

0.00701 

76 

e*  =80.1 

s,0 


Iterations 


k*  =  48.06 

S 


Figure  29:  Plots  of  eSj 0  (left)  and  ks  (right)  against  the  number  of  iterations  for 
simulations  in  which  the  other  five  parameters  are  fixed  at  values  that  have  10%-50% 
relative  error  from  the  true  values  as  given  in  (104),  and  relative  noise  of  varying 
levels  is  added  to  the  data. 
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Figure  30:  Nelder-Mead  for  the  estimation  of  eSio  and  ks  for  water.  The  other  five 
parameters  are  fixed  at  the  values  in  (104),  and  relative  random  noise  of  varying  levels 
is  added  to  the  data. 
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We  next  attempt  the  inverse  problem  with  simulated  data  to  which  we  add  con¬ 
stant  variance  noise.  We  would  like  to  compare  the  results  of  this  case  with  the  results 
obtained  from  the  inverse  problem  by  adding  relative  random  noise  to  the  simulated 
data.  We  use  the  notion  of  signal  to  noise  ratio  to  make  such  comparisons. 


6.2  Signal  to  noise  ratio 


In  analog  and  digital  communications,  the  signal  to  noise  ratio,  denoted  as  SNR, 
is  a  measure  of  the  signal  strength  relative  to  the  background  noise.  The  ratio  is 
usually  measured  in  decibels  (dB).  There  are  many  equivalent  ways  of  defining  the 
SNR  [ JN84] .  Since  we  assume  the  mean  of  our  noise  to  be  zero,  we  define  the  SNR 


as 


SNR 


10  log10 


Efcli  var(Ofc  —  E%)  J 


dB, 


(105) 


where  dB  denotes  decibels.  If  the  signal  strength  is  equal  to  the  variance  of  the 
background  noise,  i.e., 

M  M 

E  w  =  Evar(°*-Bi)  <106> 

k= 1  k= 1 


then  SNR  =  0,  and  the  signal  is  almost  unreadable,  as  the  noise  level  severely  com¬ 
petes  with  it.  Ideally,  the  strength  of  the  signal  is  greater  than  the  variance  of  the 
noise,  so  that  the  SNR  is  positive  which  results  in  the  signal  being  readable.  If  the 
strength  of  the  signal  is  less  than  the  variance  of  the  noise,  then  the  SNR  is  negative. 
Communications  engineers  always  strive  to  maximize  the  SNR  ratio.  The  SNR  ratio 
can  be  increased  by  providing  the  source  with  a  higher  level  of  signal  output  power 
if  necessary.  In  wireless  systems,  it  is  always  important  to  optimize  the  performance 
of  the  transmitting  and  receiving  antennas. 

Using  the  notation  presented  in  Section  5,  the  SNR  in  the  cases  of  relative  random 
noise  SNRr  and  constant  variance  noise  SNRC  are,  respectively,  given  as 


SNRr  =  10  log10 

SNRC  =  10  log10 


T.kL i  var (vE*krik)  J 


dB, 


gig 

Efcii  yar (/%) 


dB 


(107) 


To  compare  the  case  of  constant  variance  noise  with  the  case  of  relative  random  noise, 
we  proceed  as  follows. 

1.  As  noted  in  (79)  and  the  discussion  thereafter,  we  set  (3  =  ua.  Setting  the  SNR 
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Table  6:  Comparison  of  relative  ran¬ 
dom  noise  and  constant  variance  noise 


SNR  (dB) 

V 

%  RN 

j3  =  ua 

66 

0.0005 

0.1 

0.0118 

56 

0.0015 

0.3 

0.0353 

52 

0.0025 

0.5 

0.0588 

46 

0.005 

1.0 

0.1175 

36 

0.015 

3.0 

0.3526 

32 

0.025 

5.0 

0.5877 

26 

0.05 

10.0 

1.1754 

in  the  two  cases  to  be  equal  i.e.,  SNRr  =  SNRC,  implies 

M 


v8ly(v  El'rjk)  =  M  v2a2 


M  I  Z71*  |2 


a  = 


k= 1 


EfcLi  Iff 

M 


(108) 


2.  Thus  for  a  given  value  of  v,  we  have  [3  =  az/  and  the  corresponding  SNR  is 
calculated  as 


SNRr  =  SNRC  =  10  log 


Ef=i  Iff:  I2 


10  l  Ma2u 2 


=  101°gio(^ 


(109) 


Hence,  for  each  value  of  v  there  is  a  corresponding  value  of  /3  for  which  the  SNR  for 
relative  random  noise  is  equivalent  to  the  SNR  for  constant  variance  noise.  In  Table 
6  we  present  the  values  of  the  SNR,  /3,  and  v  that  correspond  to  each  other  in  the 
manner  described  above.  For  this  test  problem  a  =  23.5077. 

We  next  attempt  to  identify  the  two  parameters  eS)o  and  ks  by  adding  constant 
variance  noise  as  given  in  (79)  to  our  data  and  using  the  fixed  values  (104).  The 
results  for  the  corresponding  inverse  problem  are  presented  in  Table  7.  The  results  in 
this  case  are  quite  different  from  the  case  of  relative  noise  as  shown  in  Table  5.  We 
can  see  that  the  estimation  of  ks  becomes  worse  as  the  level  of  noise  increases.  The 
estimation  of  eSjo  on  the  other  hand  is  relatively  stable,  though  it  is  not  as  accurate 
as  the  estimate  obtained  in  the  case  where  relative  noise  is  added  to  the  data.  The 
estimates  also  become  worse  as  the  level  of  constant  variance  noise  increases.  In  Table 
7,  v  corresponds  to  the  scaling  factor  in  (78).  The  constant  variance  case  and  the 
relative  noise  case  for  the  same  value  of  v  have  the  same  signal  to  noise  ratio.  Figure 
31  plots  the  variation  of  eS;o  and  ks  over  all  iterations.  As  can  be  seen  here  the  value 
of  ks  increases,  away  from  the  true  value,  as  the  noise  level  increases.  Figure  32  plots 
the  details  of  the  corresponding  Nelder-Mead  Simulation. 
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Figure  31:  Plots  of  eSj 0  (left)  and  ks  (right)  against  the  number  of  iterations  for 
simulations  in  which  the  other  five  parameters  are  fixed  at  values  that  have  10%- 
50%  relative  error  from  the  true  values  as  given  in  (104),  and  constant  variance  noise 
corresponding  to  varying  levels  of  relative  noise  is  added  to  the  data. 


Table  7:  Parameter  estimation  of  eSio  and  ks  for  0%-10%  constant  variance  noise.  The 
other  five  parameters  are  fixed  at  the  values  given  in  (104) 


%  RN 

Iter 

es,0 

ks 

|J*+i  -  J£l 

J* 

l|VJA'||zx 

FC 

0.0 

38 

80.2165 

55.2945 

1.5487x  10-9 

0.4369 

0.0013 

73 

0.1 

38 

80.2114 

55.4112 

8.1590  x  10“9 

0.6973 

0.01237 

74 

0.3 

41 

80.2012 

55.6422 

9.6609  x  10“9 

2.8205 

0.001064 

77 

0.5 

40 

80.1910 

55.8728 

8.1572  x  10~9 

7.0804 

0.0005887 

77 

1.0 

41 

80.1654 

56.4531 

3.1244  x  10~9 

27.0778 

0.007561 

82 

3.0 

40 

80.0636 

58.7953 

9.7228  x  10~9 

240.6053 

0.01135 

77 

5.0 

39 

79.9621 

61.1779 

7.8005  x  10~9 

667.7927 

0.0111 

79 

10.0 

44 

79.7084 

67.3235 

3.5288  x  10~10 

2670.5115 

0.007417 

85 
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Table  8:  Confidence  Intervals  for  the  parameter  estimation  of  eS)o  and  ks  for  constant 
variance  noise  with  the  Nelder-Mead  algorithm. 


RN  =  0.0% 

es,0 

fts 

(8.02165  ±  1.8804  x  10“3)  x  10 
(5.52945  ±  1.6448  x  10"2)  x  10 

RN  =  0.1% 

es,0 

ks 

(8.02114  ±2.3718  x  10“3)  x  10 
(5.54112  ±  2.0745  x  10"2)  x  10 

RN  =  0.3% 

±s,  0 

ks 

(8.02012  ±  4.7554  x  10“3)  x  10 
(5.56422  ±4.1583  x  10"2)  x  10 

RN  =  0.5% 

£s,o 

ks 

(8.01910  ±  7.5113  x  10“3)  x  10 
(5.58728  ±  6.5668  x  10"2)  x  10 

RN  =  1.0% 

Cs,  0 

ks 

(8.01654  ±  1.4577  x  10“2)  x  10 
(5.64531  ±  1.2738  x  10’1)  x  10 

RN  =  3.0% 

±s,  0 

ks 

(8.00636  ±4.2160  x  10“2)  x  10 
(5.87953  ±  3.6783  x  lO’1)  x  10 

RN  =  5.0% 

±s,  0 

ks 

(7.99621  ±6.8206  x  10“2)  x  10 
(6.11779  ±  5.9447  x  10'1)  x  10 

RN  =  10.0% 

es,0 

ks 

(7.97084  ±  1.2723  x  lO^1)  x  10 
(6.73235  ±  1.1076)  x  10 

We  now  calculate  confidence  intervals  for  the  estimates  that  we  obtained  in  Table 
7.  The  procedure  for  calculating  these  intervals  was  outlined  in  Section  4.1.  These 
intervals  are  presented  in  Table  8  for  varying  levels  of  noise.  From  our  earlier  dis¬ 
cussion  about  confidence  intervals,  we  note  that  the  size  of  the  intervals  is  in  direct 
proportion  to  the  level  of  uncertainty  that  we  have  about  the  estimates  that  we  have 
obtained.  Thus,  the  smaller  the  confidence  interval  the  more  confident  we  are  about 
the  estimates  obtained.  As  we  can  see  from  Table  8  we  have  tighter  intervals  for  the 
estimates  of  eSi0  than  we  do  for  the  estimates  of  ks,  which  is  expected  as  our  problem 
is  more  sensitive  to  eSjo  than  it  is  to  the  pressure  coefficient  ks.  We  also  notice  that 
the  intervals  become  larger  as  the  level  of  noise  that  is  added  to  the  data  is  increased. 
These  intervals  correspond  to  a  95%  confidence  level. 
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Figure  32:  Nelder-Mead  for  the  estimation  of  eSio  and  ks  for  water.  The  other  five 
parameters  are  fixed  at  the  values  given  in  (104),  and  constant  variance  noise  corre¬ 
sponding  to  varying  levels  of  relative  noise  is  added  to  the  data. 
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Least  Squares  Objective  Function  (v  =  0.005) 


Js 


Least  Squares  Objective  Function  Cross-sections 


Figure  33:  (Top)  The  Least  Squares  Objective  Function  for  different  values  of  eSio 
and  ks.  Constant  variance  noise  with  u  =  0.005  (1%  RN)  is  added  to  the  data.  The 
other  five  parameters  were  fixed  at  their  true  values.  The  minimum  for  this  plot  is 
located  at  the  point  (80.099,49.0212),  with  a  minimum  value  of  26.712.  (Bottom) 
Cross  sections  of  the  Least  Squares  Objective  function  versus  ks  (left)  and  versus  eSi0 
(right). 


60 


The  Least  Squares  Objective  Function  (v  =  0.005) 
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Least  Squares  Function  Cross-Sections 


Figure  34:  (Top)  The  Least  Squares  Objective  Function  for  different  values  of  eSi0  and 
ks.  Constant  variance  noise  with  v  =  0.005  (1%  RN)  is  added  to  the  data.  The  other 
five  parameters  were  fixed  at  the  values  given  in  (104).  The  minimum  is  located  at  the 
point  (80.099,56.7108),  with  a  minimum  value  of  27.0837.  (Bottom)  Cross  sections 
of  the  Least  Squares  Objective  function  versus  ks  (left)  and  versus  eSio  (right). 
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Time  (sec) 

Figure  35:  Plot  of  the  Least  Squares  residual  versus  time  in  the  case  of  10%  relative 
random  noise  added  to  the  data. 


In  Figure  33  we  plot  the  least  squares  objective  function  for  values  of  eSio  that 
have  np  to  5%  relative  error  from  the  true  values,  and  values  of  ks  that  have  np  to 
20%  relative  error  from  the  true  values.  The  other  five  parameters  were  fixed  at  their 
true  values  in  the  case  of  the  top  plot  in  Figure  33  and  at  the  values  (104)  in  the  case 
of  the  bottom  plot,  respectively.  In  Figure  34  the  top  and  bottom  plots  show  the 
the  cross  sections  of  the  objective  function  against  eSjo,  and  ks  across  the  minimum 
with  respect  to  the  other  variable,  for  the  surface  plots  in  Figure  33.  Comparing  with 
Figure  25  we  observe  that  the  minimum  value  of  ks  has  moved  away  from  its  true 
value. 

In  Figures  35  we  plot  the  least  squares  residual  versus  time  in  the  case  when  10% 
relative  random  noise  is  added  to  the  data.  A  magnified  plot  of  the  least  squares 
residual  versus  time  in  the  interval  where  the  acoustic  reflection  is  observed  is  plotted 
in  Figure  36.  Finally  we  plot  the  least  squares  residual  versus  the  absolute  value 
of  the  electric  field  data  in  Figure  37  again  in  the  case  where  10%  relative  random 
noise  is  added  to  the  data.  We  notice  in  this  figure  that  the  amount  of  noise  increases 
as  the  absolute  value  of  the  electric  held  data  magnitude  increases.  This  is  in  direct 
contrast  to  the  case  when  constant  variance  noise  is  added  to  the  simulated  data.  We 
plot  analogous  plots  of  the  least  squares  residual  for  10%  added  constant  variance 
noise  in  Figures  38  and  39.  Comparing  Figures  38  and  35  we  note  the  lack  of  any 
patterns  in  the  residual  in  Figure  38  for  the  case  of  added  constant  variance  noise; 
whereas  in  Figure  35,  the  residual  follows  the  peaks  in  the  simulated  data  used  in 
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Figure  36:  A  magnified  plot  of  the  Least  Squares  residual  versus  time  in  the  interval 
where  the  acoustic  reflection  is  observed,  for  the  case  of  10%  relative  random  noise 
added  to  the  data. 


Figure  37:  Plot  of  the  Least  Squares  residual  versus  the  absolute  value  of  the  electric 
field  data  for  the  case  of  10%  relative  random  noise  added  to  the  data. 
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Figure  38:  Plot  of  the  Least  Squares  residual  versus  time  for  the  case  of  constant 
variance  noise  added  to  the  simulated  data. 


the  inverse  problem.  Again  in  Figure  37  we  note  the  fan  like  structure  of  the  residual 
with  the  residual  increasing  with  the  absoulte  value  of  the  electric  held.  In  the  case 
of  added  constant  variance  noise,  however,  we  again  note  the  lack  of  a  pattern  in 
the  residual  as  seen  in  Figure  39.  We  also  observe  that  there  are  many  points  of  the 
residual  on  the  line  Ex  =  0.  This  is  beacuse  in  the  simulated  data,  the  electric  held 
magnitude  is  almost  zero  most  of  the  time.  We  compare  the  simulated  data  with 
added  relative  random  noise  and  added  constant  variance  noise  in  Figures  40-43.  We 
notice  in  Figure  43  that  the  acoustic  reflection  is  drowned  by  the  constant  variance 
noise  with  v  =  0.05  (10%  RN),  which  is  not  the  case  with  the  relative  random  noise 
data  for  the  same  value  of  v.  Thus  the  identification  of  ks  is  more  difficult  with 
constant  variance  noise  added  to  the  simulated  data. 

We  next  attempt  the  identification  of  three  parameters,  namely  eS)o>  To  and  ks.  The 
other  four  parameters,  e^o,  cr,  «oo  and  kt  are  fixed  at  values  given  in  (104).  Tables  9 
and  10  present  the  final  estimates  and  the  details  of  the  corresponding  Nclder-Mead 
simulation.  We  observe  that  the  final  estimates  for  both  r0  and  ns  start  to  increase 
away  from  their  true  values  as  the  level  of  constant  variance  noise  is  increased.  The 
plots  of  the  least  squares  objective  function  presented  before  provide  an  explanation 
for  this.  In  Figures  45-44  we  plot  the  variation  of  the  three  parameters  over  all 
iterations  in  the  Nelder-Mead  simulation.  In  each  of  these  figures  the  dashed  line 
represents  the  true  value  of  the  corresponding  parameter.  Again,  we  observe  that 
the  estimates  for  eSio  are  stable  with  respect  to  the  level  of  constant  variance  noise, 
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Figure  39:  Plot  of  the  Least  Squares  residual  versus  the  absolute  value  of  the  electric 
held  for  the  case  of  constant  variance  noise  added  to  the  simulated  data. 


however,  the  estimates  for  tq  and  ks  increase  as  the  value  of  v  is  increased. 


Table  9:  Parameter  estimation  of  eSio  and  To  and  ks  for 
constant  variance  noise,  corresponding  to  varying  levels 
of  relative  noise  with  the  Nelder-Mead  algorithm. 


%  RN 

Iter 

es,0 

fo(xlO-“) 

ks 

0.0 

78 

80.1069 

8.1543 

48.2026 

0.1 

81 

80.1033 

8.1627 

48.3386 

0.3 

75 

80.0959 

8.1768 

48.6725 

0.5 

81 

80.0890 

8.1919 

48.9905 

1.0 

90 

80.0745 

8.2242 

49.8265 

3.0 

79 

80.0076 

8.3635 

53.1810 

5.0 

99 

79.9315 

8.5015 

56.7926 

10.0 

83 

79.7087 

8.8299 

66.2803 
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Figure  40:  Comparison  of  the  data  with  added  relative  random  noise  to  the  data  with 
constant  variance  noise  added.  The  value  of  v  in  both  the  cases  is  v  =  0.05  (10% 
RN). 


Figure  41:  Comparison  of  the  data  with  added  relative  random  noise  to  the  data  with 
constant  variance  noise  added  for  the  input  signal  at  the  center  of  the  antenna.  The 
value  of  v  in  both  the  cases  is  v  =  0.05  (10%  RN). 
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Figure  42:  Comparison  of  the  data  with  added  relative  random  noise  to  the  data  with 
constant  variance  noise  added  for  the  first  reflection  from  the  air-Debye  interface.  The 
value  of  v  in  both  the  cases  is  v  =  0.05  (10%  RN). 


Figure  43:  Comparison  of  the  data  with  added  relative  random  noise  to  the  data  with 
constant  variance  noise  added  for  the  reflection  from  the  acoustic  wave  region.  The 
value  of  v  in  both  the  cases  is  v  =  0.05  (10%  RN). 
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Figure  44:  Estimation  of  eSio  by  the  Nelder-Mead  algorithm.  We  plot  the  variation 
of  eS)o  with  the  iteration  count. 


Table  10:  Parameter  estimation  of  eS)o  and  Tq  and  ks  for 
constant  variance  noise,  corresponding  to  varying  levels 
of  relative  noise  with  the  Nelder-Mead  algorithm. 


%  RN 

||VJa||L2 

FC 

-  Jo  1 

J5 

0.0 

0.008942 

140 

8.02867xl0“9 

4.4391xl0-5 

0.1 

0.00224 

144 

4.3059  xl0“9 

0.2672 

0.3 

0.00499 

136 

5.9737  x  10-9 

2.4039 

0.5 

0.0139 

150 

3.5686x  10~9 

6.6771 

1.0 

0.001414 

161 

8.5522xl0~9 

26.7071 

3.0 

0.0277 

147 

9.0774xl0~9 

240.3553 

5.0 

0.007697 

181 

6.7874xl0~9 

667.6441 

10.0 

0.001432 

154 

6.2983xl0~9 

2670.5056 

In  the  case  of  the  estimation  of  three  parameters  we  observe  that  the  estimates 
of  ks  are  more  accurate  when  To  is  allowed  to  vary  as  opposed  to  the  case  of  two 
parameter  estimation  when  tq  is  fixed  at  a  particular  value.  Varying  To,  however, 
does  not  seem  to  make  much  difference  in  estimating  es,o-  Also,  comparing  the  in¬ 
tervals  in  Tables  8  and  11  we  find  that  we  obtain  tighter  confidence  intervals  for  the 
three  parameter  case  as  opposed  to  the  two  parameter  case,  for  the  same  level  of 
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k*  =  48.06  (-  -) 


Figure  45:  Estimation  of  ks  by  the  Nelder-Mead  algorithm.  We  plot  the  variation  of 
ks  with  the  iteration  count. 


Figure  46:  Estimation  of  r0  by  the  Nelder-Mead  algorithm.  We  plot  the  variation  of 
fo  with  the  iteration  count. 
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Table  11:  Confidence  Intervals  for  the  parameter  esti¬ 
mation  of  eS)o,To  and  ks  for  constant  variance  noise  with 
the  Nclder-Mead  algorithm. 


RN 

0.0% 

To 

ks 

(8.01069  ±  1.8738  x  HT5)  x  10 
(8.1543  ±  2.3240  x  10~4)  x  10~12 
(4.82026  ±  2.6832  x  10~4)  x  10 

RN 

0.1% 

To 

ks 

(8.01033  ±  1.4526  x  10~3)  x  10 
(8.1627  ±  1.8004  x  10~2)  x  HT12 
(4.83386  ±  2.0798  x  10~2)  x  10 

RN 

0.3% 

To 

ks 

(8.00959  ±  4.3438  x  10~3)  x  10 
(8.1768  ±  5.3824  x  10“2)  x  10”12 
(4.86725  ±  6.2220  x  HR2)  x  10 

RN 

0.5% 

To 

ks 

(8.00890  ±  7.2206  x  10“3)  x  10 
(8.1919  ±  8.9417  x  10~2)  x  10“12 
(4.89905  ±  1.0345  x  HT1)  x  10 

RN 

1.0% 

T0 

ks 

(8.00745  ±  1.4327  x  HT2)  x  10 
(8.2242  ±  1.7739  x  lO’1)  x  10“12 
(4.98265  ±  2.0555  x  HT1)  x  10 

RN 

3.0% 

T0 

ks 

(8.00076  ±  4.1809  x  10~2)  x  10 
(8.3635  ±  5.1532  x  10’1)  x  10“12 
(5.31810  ±  6.0160  x  HT1)  x  10 

RN 

5.0% 

Cs,0 

To 

ks 

(7.99315  ±  6.7761  x  10“2)  x  10 
(8.5015  ±  8.3067  x  10"1)  x  10“12 
(5.67926  ±  9.7751  x  HT1)  x  10 

RN 

10.0% 

es,o 

T0 

ks 

(7.97087  ±  1.2738  x  HT1)  x  10 
(8.8299  ±  1.5286)  x  10~12 
(6.62803  ±  1.8366)  x  10 
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added  noise.  The  reason  for  this  again  is  most  likely  the  dependence  of  the  acous¬ 
tic  reflections  on  the  value  of  tq.  Thus,  the  algorithm  for  the  inverse  problem  can 
choose  To  appropriately  so  that  the  value  of  ks  is  minimized.  As  before,  estimates  for 
all  three  parameters  become  worse  as  the  level  of  noise  added  is  increased  and  the 
corresponding  confidence  intervals  become  larger  as  seen  in  Table  11. 

In  Figure  47  we  plot  the  details  of  the  Nelder-Mead  simulation  for  the  identification 
of  the  three  parameters  eSjo,  To  and  ks.  In  Figure  48  we  plot  the  least  squares  objective 
function  for  values  of  r0  and  ks  that  differ  from  the  true  values  by  up  to  10%,  in 
the  case  of  added  constant  variance  noise  with  v  =  0.05  (10%  RN).  The  other  five 
parameters  are  fixed  at  the  values  given  in  (104)  with  eSjo  fixed  at  the  value  73.1.  We 
observe  that  the  minimum  has  moved  away  from  the  true  values.  The  fact  that  the 
acoustic  reflection  is  overshadowed  by  the  constant  variance  noise  as  shown  in  Figure 
43  explains  this  behaviour  of  the  objective  function. 

6.3  Simulation  Results:  The  Levenberg-Marquardt  method 

We  now  repeat  our  inverse  problem  with  a  different  least  squares  optimization  tech¬ 
nique.  The  Nelder  Mead  method  presented  in  the  previous  section  is  a  gradient  free 
method  based  on  the  use  of  simplices.  We  will  now  use  a  gradient  based  method, 
namely  the  Levenberg-Marquardt  method  to  solve  our  inverse  problem.  For  the 
overdetermined  least  squares  objective  function 

1  M  1 

j5(q)  =  2  \E °>  A5  q)  -  °l\2  =  2i?(q)TjR(q)’  (n°) 

Z  k=  1 

Levenberg-Marquardt,  which  is  a  trust  region  method,  adds  a  regularization  param¬ 
eter  vlm  >  0  to  the  approximate  Hessian  of  the  Gauss-Newton  method  to  obtain  a 
new  estimate  qc  =  qc  +  s,  where 

3  =  -  {ulmI  +  R'iqcfR'^y1  R\qc)TR(qc),  (HI) 

with  qc,  the  current  estimate  and  I  the  l  x  l  identity  matrix,  where  l  is  the  number  of 
parameters  that  we  are  attempting  to  estimate.  The  matrix  (vLMI  +  i?'(qc)Ti?/(qc)) 
is  positive  definite.  The  parameter  vj_,m  is  called  the  Levenberg-Marquardt  parameter. 
For  details  of  this  implementation  of  the  Levenberg-Marquardt  method  we  refer  the 
reader  to  [Kel99] . 

As  with  Nelder-Mead,  we  first  attempt  the  identification  of  all  the  seven  parame¬ 
ters,  then  we  identify  eSio  and  ks  and  finally  we  attempt  the  identification  of  eSjo,  t0 
and  ks.  Figures  49-50  plot  the  final  estimates  for  all  the  seven  parameters  and  the 
details  of  the  Levenberg-Marquardt  simulation.  The  dashed  lines  here  represent  the 
true  values  of  the  corresponding  parameter.  Tables  12  and  13  present  the  results  for 
the  inverse  problem  of  identifying  the  seven  parameters.  We  observe  how  close  the 
final  estimates  of  the  parameters  here  are  with  those  computed  by  Nelder-Mead.  We 
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Figure  47:  Nelder-Mead  for  the  three 
of  constant  variance  noise. 
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estimation  problem  for  varying  levels 
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Least  Squares  Function  Cross-Sections 


Figure  48:  (Top)  The  Least  Squares  Objective  Function  for  different  values  of  r0  and 
ks.  constant  variance  noise  with  v  =  0.005  (1%  RN)  is  added  to  the  data.  The  other 
five  parameters  were  fixed  at  the  values  given  in  (103)  with  eSjo  fixed  at  the  value  73.1. 
The  minimum  is  located  at  the  point  (8.91  x  10-12, 43.254),  with  a  minimum  value  of 
85.796.  (Bottom)  Plot  of  the  Least  Squares  Objective  function  versus  ks  (left)  and 
versus  r0  (right). 
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Norm  of  Gradient 


Least  Squares  Function  Value 


Figure  49:  Levenberg-Marquardt  for  the  seven  parameter  estimation  problem 


Figure  50:  Estimation  of  eS)o,  To  and  ks  by  the  Levenberg-Marquardt  algorithm.  We 
plot  the  variation  of  the  three  parameters  with  the  iteration  count. 
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Iterations  Iterations 


Figure  51:  Levenberg-Marquardt  for  the  seven  parameter  estimation  problem.  We 
plot  the  variation  of  the  parameters  with  the  iteration  count.  The  dashed  lines  in¬ 
dicate  the  true  values  of  the  parameters,  while  the  solid  lines  are  the  estimates  over 
the  iterations 
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Table  12:  Final  estimates  for  all  seven  parameters  for  water  using  Levenberg- 
Marquardt.  The  initial  parameter  set  has  values  that  are  5%  smaller  than  the  true 
values. 


Parameter 

True  Values 

Final  Estimates 

80.1 

80.0917 

£oo,0 

5.5 

5.2039 

T0 

8.1e  -  12 

8.0666  x  10-12 

a 

1.0  x  10"5 

9.5  x  10“6 

ks 

48.06 

48.0028 

^oo 

4.4 

4.1786 

Kr 

4.05  x  10"13 

3.8462  x  IQ'13 

Table  13:  Status  of  the  Levenberg- 
Marquardt  Simulation  at  the  final  it¬ 
eration. 


Least  squares  value 

1.0584  x  10“5 

L2  norm  of  the  gradient 

5.7248  x  10~4 

Number  of  forward  solves 

31 

Number  of  gradient  evaluations 

16 

Number  of  hessian  evaluations 

30 

can  again  make  similar  conclusions  as  we  did  in  the  case  of  the  Nelder-Mead  results. 
Only  eS)0,  t0  and  ks  seem  to  be  identifiable  with  some  accuracy,  whereas  the  problem 
appears  to  be  insensitive  to  changes  in  the  other  four  parameters.  We  next  attempt 
the  identification  of  eSio  and  ks  with  simulated  data  in  which  constant  variance  noise 
is  added.  The  values  of  the  remaining  five  parameters  are  fixed  at  those  in  (104). 
Tables  14  and  15  present  the  final  estimates  and  the  details  of  the  corresponding  LM 
simulation.  In  Figures  52  and  53  we  plot  the  details  of  the  LM  run  as  well  as  the 
variation  of  eS;o  and  ks  over  all  iterations.  The  confidence  intervals  for  these  estimates 
are  presented  in  Table  16.  Again,  we  can  make  similar  remarks  and  observations  as 
in  the  case  of  the  results  of  the  Nelder-Mead  method. 

Finally  we  attempt  to  identify  the  three  parameters  eS)o,  To  and  k,s.  Figures  55  and 
54  plot  the  details  of  the  LM  simulation  and  the  variation  of  the  parameters  over  all 
iterations.  Tables  17  and  18  present  the  final  estimates  and  details  of  the  algorithm. 
Again,  as  was  observed  in  the  results  of  the  Nelder-Mead  runs,  the  estimates  for  ks 
are  better  when  the  value  of  To  is  kept  variable  as  opposed  to  fixing  the  value  of  r0. 
Thus,  the  inverse  problem  captures  ks  more  accurately  if  the  value  of  r0  is  allowed  to 
change.  In  Table  19  we  present  the  confidence  intervals  for  the  three  parameter 
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Table  14:  Parameter  estimation  of  eS)o  and  for  constant  variance  noise,  correspond 
ing  to  varying  levels  of  relative  noise  with  the  Levenberg-Marquardt  algorithm. 


%  RN 

Iter 

£s,0 

ks 

P 

||VP||i2 

0.0 

18 

80.2166 

55.2951 

0.4369 

0.0006526 

0.5 

18 

80.191 

55.8726 

7.0804 

0.0071 

1.0 

18 

80.1654 

56.4524 

27.0778 

0.00076 

3.0 

18 

80.0636 

58.7948 

240.6053 

0.001068 

5.0 

18 

79.9621 

61.1765 

667.7927 

0.001385 

10.0 

19 

79.7084 

67.3217 

2670.5115 

0.001417 

Table  15:  Parameter  estimation  of  eS)o  and  ks  for  constant  variance  noise,  correspond 
ing  to  varying  levels  of  relative  noise  with  the  Levenberg-Marquardt  algorithm. 


%  RN 

Iter 

,P 

IIVPIU* 

F 

G 

H 

0.0 

18 

0.4369 

0.0006526 

37 

19 

36 

0.5 

18 

7.0804 

0.0071 

37 

19 

36 

1.0 

18 

27.0778 

0.00076 

37 

19 

36 

3.0 

18 

240.6053 

0.001068 

37 

19 

36 

5.0 

18 

667.7927 

0.001385 

37 

19 

36 

10.0 

19 

2670.5115 

0.001417 

39 

20 

38 

e*  =  80.1  K*  =  48.06 

s,0  s 


Figure  52:  Levenberg-Marquardt  for  the  two  parameter  estimation  problem.  Varia 
tion  of  iSfl  (left)  and  ks  (right)  over  the  iteration  count. 
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Table  16:  Confidence  Intervals  for  the  parameter  estimation  of  es  0  and  ks  for  constant 
variance  noise  with  the  Levenberg-Marquardt  algorithm. 


RN  =  0.0% 

Ks 

(8.02166  ±  1.8803  x  10“3)  x  10 
(5.52951  ±  1.6448  x  KT2)  x  10 

RN  =  0.5% 

es,0 

Ks 

(8.01910  ±  7.5113  x  10“3)  x  10 
(5.58726  ±  6.5668  x  HT2)  x  10 

RN  =  1.0% 

es,0 

Ks 

(8.01654  ±  1.4577  x  10“2)  x  10 
(5.64524  ±  1.2738  x  10"1)  x  10 

RN  =  3.0% 

Cs,  0 
Ks 

(8.00636  ±4.2160  x  10“2)  x  10 
(5.87948  ±  3.6783  x  KT1)  x  10 

RN  =  5.0% 

Cs,  0 
Ks 

(7.99621  ±6.8206  x  10“2)  x  10 
(6.11765  ±  5.9447  x  lO’1)  x  10 

RN  =  10.0% 

es,0 

Ks 

(7.97084  ±  1.2723  x  10"1)  x  10 
(6.73217  ±  1.1076)  x  10 

Norm  of  Gradient  Least  Squares  Function  Value 


Figure  53:  Estimation  of  eSio  and  ks  by  the  Levenberg-Marquardt  algorithm  for  vary¬ 
ing  levels  of  constant  variance  noise.  Details  of  the  LM  simulation. 
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Iterations  Iterations  Iterations 

Figure  54:  Estimation  of  eS)o,  tq  and  ks  by  the  Levenberg-Marquardt  algorithm  with 
varying  levels  of  added  constant  variance  noise.  We  plot  the  variation  of  the  three 
parameters  with  the  iteration  count. 
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0  5  10  15  20 

Iterations 


Least  Squares  Function  Value 
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Figure  55:  Levenberg-Marquardt  for  the  three  parameter  estimation  problem  with 
varying  levels  of  added  constant  variance  noise.  Details  of  the  LM  simulation. 
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Table  17:  Parameter  estimation  of  eSjo  and  To  and  ks  for  constant  variance  noise, 
corresponding  to  varying  levels  of  relative  noise  with  the  Levenberg-Marquardt  algo¬ 
rithm. 


%  RN 

Iter 

G,0 

To(xl0-“) 

Ks 

J5' 

0.0 

17 

80.1083 

8.1562 

48.2113 

4.4383xl0-5 

0.0009637 

1.0 

17 

80.0765 

8.2247 

49.8060 

26.7071 

0.001371 

3.0 

18 

80.0079 

8.3618 

53.1392 

240.3553 

0.001345 

5.0 

19 

79.9320 

8.4983 

56.6654 

667.6441 

0.0008044 

10.0 

20 

79.7085 

8.8330 

66.3476 

2670.5056 

0.001268 

Table  18:  Parameter  estimation  of  es  0  and  To  and 
ks  for  constant  variance  noise,  corresponding  to 
varying  levels  of  relative  noise  with  the  Levenberg- 
Marquardt  algorithm. 


%  RN 

Iter 

J5' 

livp'lb 

F 

G 

H 

0.0 

17 

4.4383  xl0“5 

0.0009637 

35 

18 

34 

1.0 

17 

26.7071 

0.001371 

35 

18 

34 

3.0 

19 

240.3553 

0.001345 

37 

19 

36 

5.0 

19 

667.6441 

0.0008044 

39 

20 

38 

10.0 

20 

2670.5056 

0.001268 

41 

21 

40 
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Table  19:  Confidence  Intervals  for  the  parameter  esti¬ 
mation  of  eS;o,To  and  ks  for  constant  variance  noise  with 
the  Levenberg-Marquardt  algorithm. 


RN 

0.0% 

To 

ks 

(8.01083  ±  1.8741  x  10“5)  x  10 
(8.1562  ±  2.3234  x  10~4)  x  10“12 
(4.82113  ±  2.6830  x  10~4)  x  10 

RN 

1.0% 

<%o 

To 

ks 

(8.00765  ±  1.4334  x  KT2)  x  10 
(8.2247  ±  1.7742  x  10'1)  x  KT12 
(4.98060  ±  2.0560  x  10”1)  x  10 

RN 

3.0% 

To 

ks 

(8.00079  ±  4.1822  x  10~2)  x  10 
(8.3618  ±  5.1552  x  10”1)  x  10”12 
(5.31392  ±6.0178  x  KT1)  x  10 

RN 

5.0% 

es,0 

To 

ks 

(7.99320  ±  6.7841  x  10~2)  x  10 
(8.4983  ±  8.3159  x  10”1)  x  10~12 
(5.66654  ±  9.7843  x  KT1)  x  10 

RN 

10.0% 

es,0 

To 

ks 

(7.97085  ±  1.2736  x  KT1)  x  10 
(8.8330  ±  1.5277)  x  10"12 
(6.63476  ±  1.8360)  x  10 

estimation  problem. 

7  A  second  test  case 

In  this  section  we  consider  the  parameter  estimation  for  a  Debye  medium  with  a  larger 
relaxation  time  than  the  value  that  was  considered  in  the  first  test  problem.  We  would 
like  to  study  the  effect  of  the  relaxation  time  r  and  the  interrogation  frequency  u>c  on 
the  ability  of  the  electromagnetic  pulse  to  penetrate  the  dielectric  material.  As  shown 
in  [BBLOO]  a  very  small  relaxation  time,  such  as  that  for  water  Tq  =  8.1  x  10~12,  will 
make  the  material  appear  comparatively  “hard”  in  the  sense  that  it  will  allow  much 
less  of  the  signal  to  penetrate  the  air- Debye  interface  at  z  —  Z\  for  a  sufficiently  high 
interrogation  frequency  like  u)c  =  n  x  1010  rad/sec.  This  small  transmitted  wave  can 
only  generate  very  little,  if  any,  reflection  at  the  second  interface  as  it  enters  the 
region  containing  the  pressure  wave.  This  demonstrates  the  importance  of  the  choice 
of  the  carrier  frequency  loc  for  effective  interrogation.  In  [BBLOO]  the  authors  have 
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demonstrated  the  strong  influence  of  the  dielectric  relaxation  time  on  an  appropriate 
choice  of  the  interrogating  frequency  for  successfully  penetrating  the  material.  In  this 
second  test  case  we  consider  a  parameter  estimation  problem  for  a  Debye  medium  with 
a  relaxation  time  Tq  =  3.16  x  1CT8  (as  used  in  [BBLOO,  ABR02]),  and  an  interrogation 
frequency  uc  =  n  x  1010  rad/sec. 

As  observed  in  Section  (4.1),  the  product  ojctq  strongly  determines  the  parameters 
in  the  problem  that  can  be  identified  accurately.  We  will  observe  that  by  changing 
the  values  of  To  and  ujc  in  this  second  test  problem,  the  parameters  that  dominate  in 
the  term  e*  in  (74)  are  no  longer  the  same  as  those  of  the  first  test  problem,  i.e.,  e*0 
and  k*.  As  will  be  demonstrated  in  Section  (7.2),  the  parameters  whose  magnitude 
dominates  the  quantity  e*  are  e^0  and  the  pressure  coefficient  k^. 

The  data  that  we  will  attempt  to  fit  is  generated  with  a  Debye  model  which  is 
defined  by  the  parameter  values 


f* 

fcs,0 

=  78.2 

f* 

too,0 

=  5.5 

T* 

To 

=  3.16  x  10 

a* 

=  1  x  10~5 

k*  =  0.05tq  =  1.58  x  10~9 
k*  =  0.6e*  0  =  46.92 
=  0.3e^0  =  1.65. 


(112) 


In  this  model,  the  value  of  Tq  is  much  higher  than  in  the  case  of  the  first  test  problem. 
As  will  be  seen  in  the  sequel,  this  value  of  Tq  leads  to  a  parameter  estimation  problem 
with  very  different  properties. 


7.1  The  forward  problem:  Simulation  results 


We  perform  similar  numerical  simulations,  as  were  done  for  the  first  test  problem, 
to  collect  simulated  data  at  the  center  of  the  antenna  and  to  demonstrate  the  effect 
of  the  acoustic  pressure  on  the  electric  field.  We  would  like  to  know  how  and  to 
what  extent  the  acoustic  pressure  can  change  the  reflected  electric  wave  from  the 
interface  at  z  =  zi-  In  Figure  56  we  plot  the  electromagnetic  source  that  is  used  in 
our  simulations.  The  form  for  the  source  is 


J s(t,  X,  z)  =  I(xl,X2)S(z)  sin (uc(t  -  3 10))  exp 

1 


in  — 


t  —  3 10 
to 

9 


(113) 


IT  X  10® 


;  u)c  =  7T  x  101U  rad/sec;  fc  =  5.0  x  10J  Hz. 


where  I(X1,X2)  is  the  Indicator  function  on  (aq,^)-  The  Fourier  spectrum  of  this  pulse 
has  even  symmetry  about  5.0  GHz.  In  Figure  57  we  plot  the  power  spectral  density 
of  the  source  defined  in  (113).  We  see  that  the  power  spectral  density  is  symmetric 
about  5.0  GHz. 

The  computational  domain  is  defined  as  follows.  We  take  A"0  =  (0,0.12),  Za  = 
(0,  0.33)  and  Z d  =  (0.33,  0.5).  The  number  of  nodes  along  the  2- axis  is  taken  to  be  450 
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and  the  number  of  nodes  along  the  x-axis  is  taken  to  be  108.  The  spatial  step  size  in 
both  the  x  and  z  directions  is  Ax  =  A z  =  h  =  0.06/54.  From  the  CFL  condition  (55) 
with  the  Courant  number  r/cN  =  1/2  we  obtain  the  time  increment  to  be  At  «  1.852 
pico  seconds.  The  central  frequency  of  the  input  source  as  described  in  (54)  is  5.0  GHz 
and  thus  the  wavelength  is  0.06  meters.  The  antenna  is  half  a  wavelength  long  and  is 
placed  at  (aq,  X2)  x  zc,  with  =  0,  aq  —  0.03  and  X2  =  0.09.  We  use  PML  layers  that 
are  half  a  wavelength  thick  on  all  fours  sides  of  the  computational  domain  as  shown 
in  Figure  3.  The  reflections  of  the  electromagnetic  pulse  at  the  air-Debye  interface 
and  from  the  acoustic  pressure  wave  are  recorded  at  the  center  of  the  antenna  (xc,  zc), 
with  xc  =  0.06,  at  every  time  step.  This  data  will  be  used  as  observations  for  the 
parameter  identification  problem  to  be  presented  later.  The  component  of  the  electric 
held  that  is  of  interest  here  is  the  Ex  component.  Thus  our  data  is  the  set 


E(q*) 

q* 


{Ex{t  =  nAt,  xc,  zc]  q*)};/=1 

/  *  *  *  *  *  *  * 
Vs,0>  eoo,0>  T0  1®  1  Koo>  ^ r ) 


(114) 


The  windowed  acoustic  pressure  wave  as  defined  in  (49)  has  the  parameter  values 
tup  =  4.07T  x  104  rad/sec,  cp  =  667.67  m/s,  and  thus,  \p  =  0.033.  The  location  of  the 
pressure  region  is  in  the  interval  (^2,  £2  +  Ap)  =  (0.44,  0.47). 

In  Figures  58  -61  we  plot  the  Ex  held  magnitude  in  the  plane  containing  the  center 
of  the  antenna  versus  depth  along  the  z  axis.  Figure  58  depicts  the  electromagnetic 
wave  penetrating  the  Debye  medium.  In  59  we  see  the  reflection  of  the  electromagnetic 
wave  from  the  Debye  medium  moving  towards  the  antenna  and  the  Brillouin  precursor 
propagating  in  the  Debye  medium.  In  Figure  60  we  observe  the  reflection  from  the 
region  containing  the  acoustic  pressure  moving  towards  the  air-Debye  interface  and 
the  transmitted  part  of  the  electromagnetic  source  travelling  into  the  absorbing  layer. 
The  reflection  from  the  acoustic  pressure  crosses  the  air-Debye  interface  in  Figure  61. 
Here  we  see  another  reflected  wave  from  the  air-Debye  interface  that  moves  towards 
the  right.  In  Figures  62  and  63  we  plot  the  Ex  held  magnitude  recorded  at  the  center 
of  the  antenna  versus  time,  which  shows  the  electromagnetic  source,  the  reflection 
off  the  air-Debye  interface  and  the  reflection  from  the  region  containing  the  acoustic 
pressure  wave.  As  seen  here  the  magnitude  of  the  reflection  from  the  acoustic  pressure 
is  many  orders  of  magnitude  smaller  that  the  initial  electromagnetic  source  as  well 
the  reflection  from  the  air-Debye  interface. 


7.2  Sensitivity  analysis 

Here  we  examine  the  system  dynamics  as  the  parameters  vary.  We  are  interested  in 
the  changes  produced  by  the  coefficents  of  pressure  in  the  polarization,  namely  the 
parameters  ks,  and  kt.  In  Figure  65  we  plot  the  difference  in  the  Ex  component 
of  the  electric  held  between  the  simulation  with  parameter  values  given  in  (112) 
and  simulations  for  which  the  value  of  k0 0  is  0,  0. and  ClAe^o.  As  expected 
the  difference  is  seen  in  the  acoustic  reflection  that  is  observed  at  the  center  of  the 
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Figure  56:  Forward  Simulation  for  a  Debye  medium  with  parameters  given  in  (112). 
The  input  electromagnetic  source  as  defined  in  (113). 


antenna.  In  Figure  66  we  plot  the  difference  in  the  Ex  component  of  the  electric 
held  between  the  simulation  with  parameter  values  given  in  (112)  and  simulations  for 
which  the  value  of  ks  is  0,  and  0.3eSjo-  In  Figure  67  we  plot  the  difference  in  the  Ex 
component  of  the  electric  held  between  the  simulation  with  parameter  values  given  in 
(112)  and  simulations  for  which  the  value  of  nT  is  0  and  0.1eSio-  Again,  in  both  cases 
the  difference  is  seen  in  the  acoustic  reflection  that  is  observed  at  the  center  of  the 
antenna.  We  observe  that  the  difference  in  the  Ex  held  in  the  case  of  Figures  66  and 
67  is  one  or  two  orders  of  magnitude  smaller  than  that  in  Figure  65.  This  indicates 
that  the  coefhcent  that  is  most  signihcant  for  this  test  problem  and  which  will  be 
estimated  more  accurately  than  the  others  is  the  pressure  coefficient  in  e^.  In 
the  hrst  test  problem  ks,  the  pressure  coefficient  in  es  was  the  most  signihcant.  So 
we  see  a  very  different  behaviour  of  the  acoustic  reflection  in  this  case. 

This  observation,  as  was  done  for  the  hrst  test  problem,  can  also  be  deduced  from 
an  analysis  of  (40).  In  the  simulations  for  this  test  problem  the  outgoing  and  rehected 
radiation  will  be  dominated  by  frequencies  near  the  center  frequency  5.0  GHz.  Thus, 
e*  will  be  dominated  by  frequencies  near  5.0  GHz.  In  this  problem  cuTq  ~  0(  102), 
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Figure  57:  Fourier  transform  of  the  source  defined  in  (113).  The  transform  is  centered 
around  the  central  frequency  of  to  =  n  x  1010  rad/sec. 


and  with  e0  =  8  x  10  12  we  have  (compare  with  the  estimates  of  (75)) 


ss,0 


1  +  j^T0* 

* 

1+j utJ  °°’° 

a* 
j^e0 


<3(10-2<0) 

fa  0{ KT1 

0(C,o) 

fa  0(10), 

£>(10V) 

fa  o{ KT3 

(115) 


Thus  we  can  see  that  e*  will  be  most  sensitive  to  the  infinite  frequency  permittivity 
Coc  o  and  the  effects  of  es>0  and  r0  and  o  will  not  be  as  pronounced.  Also,  this  implies 
that  e*  will  be  more  sensitive  to  the  pressure  coefficient  than  to  the  coefficients 
ks  and  kt,  which  is  the  same  observation  that  was  made  from  Figures  65,  66  and  67 
We  next  study  the  effect  that  the  acoustic  speed  and  frequency  have  on  the  am¬ 
plitude  of  the  reflection  from  the  acoustic  pressure  region  as  was  done  in  the  first 
test  problem.  In  Figure  68  we  plot  the  acoustic  reflection  observed  at  the  center  of 
the  antenna  for  different  values  of  the  acoustic  frequency  and  in  Figure  69  we  plot 
the  acoustic  reflection  for  different  values  of  the  acoustic  speed.  As  noted  before,  by 
changing  the  speed  or  the  frequency,  the  wavelength  Xp  changes  and  thus  the  size 
of  the  interval  (z2,Z2  +  Ap)  in  which  the  pressure  wave  is  generated.  The  amplitude 
of  the  acoustic  reflection  appears  to  increase  as  the  acoustic  frequency  is  increased 
and  decreases  as  the  acoustic  speed  increases,  the  opposite  behaviour  to  what  was 
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t  =  1 .85  x  1 0'9  sec 


Figure  58:  Forward  Simulation  for  a  Debye  medium  with  parameters  given  in  (112). 
The  Ex  field  magnitude  along  the  plane  containing  the  center  of  the  antenna  versus 
depth.  We  see  here  the  electromagnetic  source  penetrating  the  Debye  medium. 


t  =  2.78  x  10-9  sec 


Figure  59:  Forward  Simulation  for  a  Debye  medium  with  parameters  given  in  (112). 
The  Ex  field  magnitude  along  the  plane  containing  the  center  of  the  antenna  versus 
depth.  The  plot  shows  the  reflection  of  the  electromagnetic  source  off  the  air-Debye 
interface  travelling  towards  the  antenna.  We  also  see  the  brillouin  precursor  propa¬ 
gating  inside  the  Debye  medium. 
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Figure  60:  Forward  Simulation  for  a  Debye  medium  with  parameters  given  in  (112). 
The  Ex  field  magnitude  along  the  plane  containing  the  center  of  the  antenna  versus 
depth.  Here  we  see  the  reflection  from  the  acoustic  pressure  wave  moving  towards 
the  left. 


t  =  4.63  x  10“9  sec 


z  (meters) 


Figure  61:  Forward  Simulation  for  a  Debye  medium  with  parameters  given  in  (112). 
The  Ex  field  magnitude  along  the  plane  containing  the  center  of  the  antenna.  The 
reflection  from  the  acoustic  pressure  wave  crosses  the  air-Debye  interface.  A  secondary 
reflection  off  this  interface  travels  towards  the  right. 
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Figure  62:  Forward  Simulation  for  a  Debye  medium  with  parameters  given  in  (112). 
The  Ex  field  magnitude  versus  time.  We  see  three  sets  of  observations.  The  electro¬ 
magnetic  source,  the  reflection  off  the  air-Debye  interface  and  the  rclection  from  the 
acoustic  pressure  wave. 
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Figure  63:  Forward  Simulation  for  a  Debye  medium  with  parameters  given  in  (112). 
The  Ex  held  magnitude  versus  time.  Here  we  magnify  the  two  sets  of  reflections,  one 
from  the  air-Debye  interface  and  from  the  acoustic  pressure  wave. 


Figure  64:  Fourier  transform  of  the  electric  field  data  in  Figure  62.  The  transform  is 
centered  around  the  central  frequency  of  a j  —  107r  x  109. 


Figure  65:  Forward  Simulation  for  a  Debye  medium  with  parameters  given  in  (112). 
The  Ex  field  magnitude  measured  at  the  center  of  the  antenna  for  different  values  of 
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Figure  66:  Difference  in  the  Ex  field  magnitude  measured  at  the  center  of  the  antenna 
between  the  simulation  for  parameters  given  in  (112)  and  simulations  with  different 
values  of  ks. 


observed  in  the  first  test  problem.  We  note  here  again  that  as  onr  windowed  pressure 
wave  contains  only  one  wavelength  of  the  sinusoid  it  is  difficult  to  predict  how  the 
reflections  will  behave  as  the  speed  or  the  frequency  is  changed. 

7.3  Simulation  results  for  the  inverse  problem:  The  Leven- 
berg  Marquardt  method 

In  this  section  we  present  results  for  the  inverse  problem  of  parameter  estimation  for 
the  second  test  problem.  We  do  this  first  with  the  Levenberg-Marquardt  algorithm 
and  later  with  the  Nelder-Mead  method.  We  first  attempt  to  estimate  the  three 
parameters  e^o,  To  and  the  pressure  coefficient  Kqo,  of  e^.  The  values  of  the  remaining 
parameters  are  fixed  at 

e*0  =71.09  (relative  static  permittivity), 

eoo  o  =5.5  (relative  high  frequency  permittivity), 

Tq  =  3.48  x  10-8  seconds,  (H6) 

cr*  =  1.5  x  10“5  (mhos/meter), 
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Figure  67:  Difference  in  the  Ex  field  magnitude  measured  at  the  center  of  the  antenna 
between  the  simulation  for  parameters  given  in  (112)  and  simulations  with  different 
values  of  nT. 


< 

=  42.654  (pressure  coefficient  in  e*), 

=  1.65  (pressure  coefficient  in  e^), 

(117) 

< 

=  1.74  x  10~9  (pressure  coefficient  in  r*). 

The  results  are  tabulated  in  Tables  20  and  21  and  presented  in  Figures  70-73.  We 
have  added  constant  variance  noise  with  different  values  of  the  scale  parameter  v  to 
the  data.  We  note  from  Table  20  that  the  estimation  of  e^o  is  relatively  stable  over 
all  values  of  v  that  we  have  selected;  whereas  the  estimation  of  becomes  inaccurate 
as  the  value  of  v  is  increased.  Unlike  the  previous  test  problem,  the  estimation  of  t$  in 
this  case  is  not  good.  Table  21  presents  the  number  of  function,  gradient  and  hessian 
evaluations  that  are  required  for  the  algorithm  to  converge,  where  the  convergence 
criteria  used  here  is  1 1  VJf  1 1/|  |  VJg  1 1  <  1CT6. 

In  Table  22  we  present  the  confidence  intervals  for  the  three  parameter  estimation 
problem  with  results  tabulated  in  Table  20.  Again,  we  note  that  the  confidence 
intervals  become  larger  as  the  level  of  noise  increases.  From  Table  22  we  observe  that 
the  confidence  level  in  the  estimates  and  k ^  is  higher  than  the  confidence  level  in 
fo,  as  is  also  observed  from  the  final  estimates  for  r0. 

We  now  attempt  to  estimate  the  two  parameters  600,0  and  its  pressure  coefficient 
ftoo.  We  fix  the  value  of  the  relaxation  time  tq  to  be 

r0  =  3.48  x  10“8.  (118) 
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Figure  68:  Plot  of  the  magnitude  of  the  acoustic  reflection  against  time  for  different 
values  of  the  acoustic  frequency.  The  other  parameters  are  fixed  at  their  true  values 
given  in  (112). 
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Figure  69:  Plot  of  the  magnitude  of  the  Ex  component  of  the  electric  field  against 
time  for  different  values  of  the  acoustic  speed.  The  other  parameters  are  fixed  at  their 
true  values  given  in  (112). 
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Table  20:  Parameter  estimation  of  e^o,  To  and  for  constant  variance  noise,  corre¬ 
sponding  to  varying  levels  of  relative  noise  with  the  Levenberg-Marquardt  algorithm. 


%  RN 

Iter 

^oo,0 

T0 

^OO 

J8 

l|VJs|k* 

0.0 

32 

5.4999 

2.8522  x  10~8 

1.6496 

2.1605  x  10~6 

0.0002746 

0.1 

32 

5.4998 

2.8761  x  10~8 

1.6463 

0.2749 

0.0002723 

0.3 

32 

5.4996 

2.9253  x  10~8 

1.6395 

2.4743 

0.0003417 

0.5 

32 

5.4993 

2.9764  x  10~8 

1.6328 

6.8730 

0.0004703 

1.0 

34 

5.4985 

3.1135  x  10~8 

1.6160 

27.4921 

0.0002261 

3.0 

39 

5.4956 

3.8493  x  10~8 

1.5467 

247.4282 

0.0002441 

5.0 

45 

5.4926 

5.1264  x  10~8 

1.4753 

687.2987 

0.0003073 

10.0 

57 

5.4852 

4.1950  x  10~7 

1.2957 

2749.1770 

0.00007404 

Table  21:  Parameter  estimation  of  e^o,  t0  and  for  constant  variance  noise,  corre¬ 
sponding  to  varying  levels  of  relative  noise  with  the  Levenberg-Marquardt  algorithm. 


%  RN 

Iter 

J8 

||vj8||L2 

F 

G 

H 

0.0 

32 

2.1605  x  10~6 

0.0002746 

71 

33 

70 

0.1 

32 

0.2749 

0.0002723 

71 

33 

70 

0.3 

32 

2.4743 

0.0003417 

71 

33 

70 

0.5 

32 

6.8730 

0.0004703 

71 

33 

70 

1.0 

34 

27.4921 

0.0002261 

76 

35 

75 

3.0 

39 

247.4282 

0.0002441 

86 

40 

85 

5.0 

32 

687.2987 

0.0003073 

98 

46 

97 

10.0 

57 

2749.1770 

0.00007404 
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Figure  70:  Parameter  estimation  of  e^o,  r0  and  k0 q  using  the  Levenberg-Marquardt 
algorithm  with  added  constant  variance  noise  for  various  values  of  v. 


Figure  71:  Levenberg-Marquardt  for  the  three  parameter  estimation  problem  with 
added  constant  variance  noise.  Variation  of  the  estimate  e^o  with  the  iteration 
count. 
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Figure  72:  Levenberg-Marquardt  for  the  three  parameter  estimation  problem  with 
added  constant  variance  noise.  We  plot  the  variation  of  the  estimate  tq  with  iteration 
count  (top)  and  a  magnified  plot  of  the  same  (bottom). 
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Table  22:  Confidence  Intervals  for  the  parameter  estima¬ 
tion  of  600,0,  To  and  /too  with  the  Levenberg-Marquardt 
algorithm  for  varying  levels  of  constant  variance  noise 
corresponding  to  0%-10%  relative  noise. 


RN  =  0.0% 

^00, 0 

T0 

^00 

5.4999  ±  9.8421  x  HT7 
(2.8522  ±  1.9957  x  10"4)  x  10“8 
1.6496  ±  2.8669  x  10"5 

RN  =  0.1% 

^00, 0 

To 

^00 

5.4998  ±  3.5104  x  10”4 
(2.8761  ±  7.1807  x  10"2)  x  10“8 
1.6463  ±  1.0238  x  10~2 

RN  =  0.3% 

^00, 0 

T0 

^00 

5.4996  ±  1.0531  x  10~3 
(2.9253  ±  2.1927  x  10”1)  x  10“8 
1.6395  ±  3.0784  x  10“2 

RN  =  0.5% 

^00, 0 

T0 

^00 

5.4993  ±  1.7548  x  10“3 
(2.9764  ±  3.7209  x  lO’1)  x  10~8 
1.6328  ±5.1417  x  10“2 

RN  =  1.0% 

^00, 0 

T0 

^00 

5.4985  ±  3.5076  x  10“3 
(3.1135  ±  7.7975  x  10”1)  x  10~8 
1.6160  ±  1.0340  x  10"1 

RN  =  3.0% 

^00, 0 

T0 

^00 

5.4956  ±  1.0504  x  10“2 
(3.8493  ±  2.9130)  x  10~8 
1.5467  ±3.1752  x  10”1 

RN  =  5.0% 

^00, 0 

T0 

^00 

5.4926  ±  1.7460  x  10“2 
(5.1264  ±6.5133)  x  10~8 
1.4753  ±  5.4261  x  10’1 

RN  =  10.0% 

^00, 0 

T0 

^00 

5.4852  ±  3.4394  x  10“2 
(4.1950  ±  10.3835)  x  10~8 
1.2957  ±  11.5944  x  10"1 
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Figure  73:  Levenberg-Marquardt  for  the  three  parameter  estimation  problem  with 
added  constant  variance  noise.  We  plot  the  variation  of  the  estimate  k ^  with  the 
iteration  count. 


In  Figures  74and  75  we  plot  the  two  parameter  estimation  inverse  problem  with  no 
added  noise  to  the  simulated  data.  In  Figures  76  and  77  we  plot  the  two  parameter 
estimation  inverse  problem  with  added  constant  variance  noise  to  the  simulated  data. 
In  Tables  23  and  24  we  present  the  final  estimates  and  details  of  the  Levenberg- 
Marquardt  simulations  with  varying  levels  of  constant  variance  noise  added  to  the 
data  that  is  to  be  fit.  In  Table  25  we  calculate  the  confidence  intervals  for  the 
two  parameter  estimation  problem  of  estimating  e^o  and  with  added  constant 
variance  noise  to  the  simulated  data.  In  Figure  78  we  plot  the  least  squares 

objective  function  for  different  values  of  e^o  and  Kqo.  There  is  no  noise  added  to  the 
data  for  these  plots.  The  other  five  parameters  are  fixed  at  their  true  values  as  given 
in  (112).  We  see  possibilities  of  local  minima  in  these  plots.  This  test  problem  has 
many  local  minima  whose  evidence  will  be  seen  clearly  in  later  plots.  In  Figure  79  we 
plot  the  least  squares  objective  function  for  various  values  of  the  two  parameters  600,0 
and  Koo.  Again  we  have  not  added  any  noise  to  the  data  for  these  plots.  The  other 
five  parameters  were  fixed  at  values  given  in  (116)  and  (117)  except  for  tq  which  was 
fixed  at  its  true  value  as  given  in  (112). 
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Figure  74:  Parameter  estimation  of  600,0  and  Koo  with  no  added  noise  using  the 
Levenberg-Marquardt  algorithm.  The  plot  shows  the  details  of  the  Levenberg- 
Marquardt  algorithm. 


eoo0=5-5  7  =  1.65 


Figure  75:  Levenberg-Marquardt  for  the  two  parameter  estimation  problem  with  no 
added  noise.  We  plot  the  variation  of  e oo,0  (left)  and  k^  (right)  versus  the  iteration 
count. 
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Figure  76:  Parameter  estimation  of  e^o  and  with  the  Levenberg-Marquardt  al¬ 
gorithm  for  varying  levels  of  constant  variance  noise  (v),  corresponding  to  0%-10% 
relative  noise.  We  plot  the  variation  of  e^o  (left)  and  koo  (right)  versus  the  iteration 
count. 


Norm  of  Gradient  Least  Squares  Function  Value 


Figure  77:  Details  of  the  Levenberg-Marquardt  simulations  for  the  two  parameter 
estimation  problem  with  varying  levels  (y)  of  constant  variance  noise  added  to  the 
simulated  data  as  tabulated  in  Table  23. 
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Table  23:  Parameter  estimation  of  600,0  and  Koo  with 
the  Levenberg-Marquardt  algorithm  for  varying  levels  of 
constant  variance  noise,  corresponding  to  0%-10%  rela¬ 
tive  noise. 


%  RN 

Iter 

^00, 0 

^OO 

R 

l|VJs||y 

0.1 

30 

5.4988 

1.5708 

0.2924 

0.0003925 

0.3 

30 

5.4986 

1.5714 

2.4885 

0.0003861 

0.5 

30 

5.4984 

1.5721 

6.8843 

0.0003768 

1.0 

29 

5.4979 

1.5738 

27.4975 

0.0007438 

3.0 

44 

5.4961 

1.5806 

247.4317 

0.0003701 

5.0 

23 

5.4943 

1.5873 

687.3371 

0.0002948 

10.0 

24 

4.9006 

-1.0116 

2800.2627 

0.0006955 

Table  24:  Parameter  estimation  of  e^o  and  with 
the  Levenberg-Marquardt  algorithm  for  varying  levels 
of  constant  variance  noise,  corresponding  to  0%-10%  of 
relative  noise. 


%  RN 

Iter 

R 

IIVR'IIls 

F 

G 

H 

0.1 

30 

0.2924 

0.0003925 

65 

31 

64 

0.3 

30 

2.4885 

0.0003861 

65 

31 

64 

0.5 

30 

6.8843 

0.0003768 

65 

31 

64 

1.0 

29 

27.4975 

0.0007438 

63 

30 

62 

3.0 

44 

247.4317 

0.0003701 

91 

45 

90 

5.0 

23 

687.3371 

0.0002948 

49 

24 

48 

10.0 

24 

2800.2627 

0.0006955 

52 

25 

51 
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Table  25:  Confidence  Intervals  for  the  pa¬ 
rameter  estimation  of  e^o  and  with  the 
Levenberg-Marquardt  algorithm  for  varying 
levels  {y)  of  constant  variance  noise  corre¬ 
sponding  to  0%-10%  relative  noise. 


RN  =  0.1% 

^oo,0 

^oo 

5.4988  ±  3.3500  x  10”4 
1.5708  ±  2.0666  x  10~3 

RN  =  0.3% 

^oo,0 

^oo 

5.4986  ±  9.7663  x  10"4 
1.5714  ±6.0237  x  10~3 

RN  =  0.5% 

^oo,0 

^oo 

5.4984  ±  1.6231  x  10“3 
1.5721  ±  1.0009  x  10-2 

RN  =  1.0% 

^oo,0 

^oo 

5.4979  ±  3.2375  x  10“3 
1.5738  ±  1.9955  x  10“2 

RN  =  3.0% 

^oo,0 

^oo 

5.4961  ±  9.6379  x  10“3 
1.5806  ±  5.9286  x  10“2 

RN  =  5.0% 

^oo,0 

^oo 

5.4943  ±  1.5943  x  10"2 
1.587  ±9.7880  x  10~2 
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600 


The  Least  Squares  Objective  Function 


Least  Squares  Objective  Function  Cross-Sections 


Figure  78:  (Top)  The  Least  Squares  Objective  Function  for  different  values  of  €00,0 
and  FCoc  with  no  noise  added  to  the  data.  The  other  five  parameters  were  fixed  at 
their  true  values  as  given  in  (112).  The  minimum  is  located  at  the  point  (5.5, 1.65), 
with  a  minimum  value  of  0.  (Bottom)  Cross  Sections  of  the  Least  Squares  Objective 
function  versus  Koo  (right)  and  versus  e^o  (left). 


103 


In  Figure  80  we  plot  the  least  squares  objective  function  for  a  range  of  values  of 
600,0  and  k oo,  again  with  no  added  noise.  In  this  case  we  have  not  added  any  noise  to 
the  data.  The  other  five  parameters  were  fixed  at  values  given  in  (116),  and  (117), 
and  To  fixed  at  the  value  given  in  (118)  that  have  10-50%  relative  error  as  compared 
with  their  true  values  in  (112). 

7.4  Simulation  results  for  the  inverse  problem:  The  Nelder- 
Mead  method 

We  first  attempt  to  estimate  the  two  parameters  e^o  and  its  pressure  coefficient 
We  add  constant  variance  noise  to  our  simulated  data.  The  other  five  parameters  are 
fixed  at  their  values  given  in  (116),  and  (117)  with  r0  fixed  at  the  value  given  in  (118). 
Figures  81  and  82  plot  the  details  of  the  Nelder-Mead  simulations  and  the  variation  of 
the  two  parameters  over  all  iterations.  Table  26  displays  the  final  estimates  and  the 
details  of  the  Nelder-Mead  algorithm  for  the  two  parameter  estimation  problem  with 
added  constant  variance  noise.  We  note  that  the  Nelder-Mead  algorithm  does  not 
converge  to  the  optimal  values  of  either  of  the  two  parameters.  Instead  it  appears  to 
converge  to  a  possible  local  minimum.  We  demonstrate  this  fact  by  plotting  the  least 
squares  objective  function  for  the  case  of  added  constant  variance  noise  corresponding 
to  v  —  0.005  (1%  RN)  in  the  Figure  83.  This  figure  clearly  demonstrates  the  presence 
of  many  local  optima.  The  reason  for  the  presence  of  local  optima  is  the  fact  that 
we  are  changing  e^o  inside  the  Debye  media.  This  changes  the  speed  of  propagation 
of  the  electromagnetic  pulse  inside  the  Debye  media  as  the  maximum  wavespeed  in 
the  dispersive  dielectric  is  Co/yT^jo  [Pet94] .  A  similar  phenomenon  was  observed  in 
[BGW03]  where  the  authors  constructed  a  modified  least  squares  objective  function 
to  avoid  the  occurence  of  local  optima.  Another  option  to  solve  this  problem  would  be 
to  resort  to  using  a  global  optimization  method  to  solve  the  inverse  problem.  We  note 
that  the  Levenberg-Marquardt  method  did  not  seem  to  converge  to  a  local  minima. 
However,  this  may  be  a  fortunate  occurence  in  which  case  we  can  always  find  some 
initial  guesses  that  would  cause  such  local  optimization  algorithms  to  converge  to  a 
local  minimum  instead  of  the  global  minimum. 

In  Table  27  we  calculate  the  confidence  intervals  for  the  two  parameter  estimation 
problem.  We  note  that  the  intervals  are  larger  as  compared  with  the  results  of  the 
Levenberg-Marquardt  method.  We  also  plot  the  least  squares  objective  function 
for  different  values  of  600,0  and  without  adding  noise  to  the  data.  Again  we  see 
the  presence  of  local  optima  as  seen  in  Figures  84.  To  see  how  the  interrogating 
frequency  relates  to  the  presence  of  local  minima  we  repeat  the  inverse  problem  by 
lowering  the  interrogating  frequency  to  3.0  GHz  from  5.0  GHz.  In  Figure  85  We 
observe  that  in  this  case  the  distance  between  successive  minima  has  increased.  Thus 
by  lowering  the  interrogating  frequency  we  can  push  the  local  minima  away  from  each 
other  and  thius  reduce  the  number  of  local  minima  that  would  occur  in  an  interval  of 
interest.  We  next  attempt  to  estimate  the  three  parameters  r0,  600,0  and  its  pressure 
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Least  Squares  Objective  Function 


5.3  5.5  5.7  5.9  1.4  1.6  1.8  2 


Figure  79:  (Top)  The  Least  Squares  Objective  Function  for  different  values  of  600,0  and 
/too  with  no  noise  added  to  the  data.  The  other  five  parameters  were  fixed  at  values 
given  in  (116),  and  (117)  except  for  tq  which  was  fixed  at  its  true  value.  The  minimum 
is  located  at  the  point  (5.5, 1.617),  with  a  minimum  value  of  0.0147.  (Bottom)  Cross 
Sections  of  the  Least  Squares  Objective  function  versus  (right)  and  versus  600,0 
(left). 
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Least  Squares  Objective  Function 

40CK 


Least  Squares  Objective  Function  Cross-sections 


Figure  80:  (Top)  The  Least  Squares  Objective  Function  for  different  values  of  €00,0 
and  Koo  with  no  noise  added  to  the  data.  The  other  five  parameters  were  fixed  at 
values  given  in  (116),  and  (117),  and  r0  fixed  at  the  value  given  in  (118)  that  have 
10-50%  relative  error  as  compared  with  their  true  values.  The  minimum  is  located  at 
the  point  (5.5, 1.584),  with  a  minimum  value  of  0.03595.  (Bottom)  Cross  Sections  of 
the  Least  Squares  Objective  function  versus  Kqo  (right)  and  versus  e^o  (left). 
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Figure  81:  Details  of  the  Nelder-Mead  algorithms  for  the  two  parameter  estimation 
of  600,0  and  Koo  for  varying  levels  (z/)  of  constant  variance  noise  added  to  the  data  as 
tabulated  in  Table  26. 
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4.8 


Figure  82:  Variation  in  the  estimates  e00j 0  (left)  and  k ^ 
(right)  over  all  iterations  in  the  Nelder-Mead  algorithm 
for  varying  levels  (z/)  of  constant  variance  noise  added 
to  the  data  as  tabulated  in  Table  26. 


Table  26:  Parameter  estimation  of  £00,0  and  k0 0  for  varying  levels  of  constant  variance 
noise  corresponding  to  0%-10%  relative  noise. 


%  RN 

Iter 

Fx),o 

^OO 

J* 

||VJy||L2 

FC 

0.1 

39 

4.5226 

1.1231 

8.5551  x  10“9 

113.3545 

0.0551 

79 

0.3 

40 

4.5224 

1.1232 

3.6791  x  10~9 

114.9105 

0.0440 

81 

0.5 

42 

4.5221 

1.1233 

1.6051  x  10~9 

118.6662 

0.0209 

85 

1.0 

42 

4.5215 

1.1236 

8.8755  x  10“9 

137.6788 

0.0610 

83 

3.0 

42 

4.5189 

1.1246 

5.7953  x  10~9 

351.2076 

0.0154 

83 

5.0 

40 

4.5164 

1.1257 

9.4709  x  10“9 

784.7016 

0.0439 

79 

10.0 

42 

4.5102 

1.1284 

3.7162  x  10~9 

2830.7864 

0.0448 

86 

108 


Table  27:  Confidence  Intervals  for  the  parameter  estima¬ 
tion  of  Cqc  g  and  with  the  Nelder-Mead  algorithm  for 
varying  levels  of  constant  variance  noise  corresponding 
to  0%-10%  relative  noise. 


RN  =  0.1% 

^oo,0 

^oo 

4.5226  ±  7.6273  x  10~3 
1.1231  ±4.7312  x  10~2 

RN  =  0.3% 

^oo,0 
^  oo 

4.5224  ±  7.6781  x  10~3 
1.1232  ±4.7624  x  10~2 

RN  =  0.5% 

^oo,0 

^oo 

4.5221  ±  7.8008  x  10~3 
1.1233  ±4.8383  x  10~2 

RN  =  1.0% 

^oo,0 

^oo 

4.5215  ±  8.3986  x  10~3 
1.1236  ±5.2073  x  10~2 

RN  =  3.0% 

^oo,0 

4.5189  ±  1.3385  x  10~2 
1.1246  ±  8.2935  x  10~2 

RN  =  5.0% 

^oo,0 

^oo 

4.5164  ±  1.9965  x  10“2 
1.1257  ±  1.2361  x  lO’1 

RN  =  10.0% 

^oo,0 

^oo 

4.5102  ±  3.7722  x  10“2 
1.1284  ±2.3306  x  10”1 
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Least  Sauares  Objective  Function  (v  =  0.005) 


Least  Squares  Objective  Function  Cross-sections 


Figure  83:  (Top)  The  Least  Squares  Objective  Function  for  different  values  of  €00,0 
and  /too  for  constant  variance  noise  with  v  =  0.005  (1%  RN)added  to  the  data.  The 
other  five  parameters  were  fixed  at  values  given  in  (116),  and  (117).  The  minimum 
is  located  at  the  point  (5.5,1.584),  with  a  minimum  value  of  27.5152.  However,  we 
note  the  presence  of  a  local  minimum  at  the  point  (4.51, 1.122)  with  a  minimum  least 
squares  value  of  138.0849.  (Bottom)  Cross  sections  of  the  Least  Squares  Objective 
function  versus  k0 0  (right)  and  versus  e^o  (left). 
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Least  Squares  Objective  Function 


Figure  84:  (Top)  The  Least  Squares  Objective  Function  for  different  values  of  600,0 
and  Koo  with  no  noise  added  to  the  data.  The  other  five  parameters  were  fixed  at 
values  given  in  (116),  and  (117),  and  To  fixed  at  the  value  in  (118).  The  minimum 
is  located  at  the  point  (5.5, 1.58),  with  a  minimum  value  of  0.036.  However,  we  note 
the  presence  of  many  local  minimum.  (Bottom)  Cross  Sections  of  the  Least  Squares 
Objective  function  versus  k0 0  (right)  and  versus  600,0  (left). 
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Least  Squares  Objective  Function:  Cross-Sections 


Figure  85:  (Top)  The  Least  Squares  Objective  Function  for  different  values  of  e^o 
and  Koo  with  no  noise  added  to  the  data.  The  other  five  parameters  were  fixed  at 
values  given  in  (116),  and  (117).  However  the  electromagnetic  input  source  had  a 
center  frequency  of  3.0  GHz  as  opposed  to  5.0  GHz.  The  minimum  is  located  at  the 
point  (5.5, 1.683),  with  a  minimum  value  of  6.91.  We  again  note  the  presence  of  local 
minimum.  (Bottom)  Cross  sections  of  the  Least  Squares  Objective  function  versus 
/too  (right)  and  versus  £00,0  (left). 


112 


Table  28:  Parameter  estimation  of  es>o,  To,  and  ks  for  0%,  1%,  3%  and  5%  constant 
variance  noise  (test  2). 


%  RN 

Iter 

7oo,0 

r0(xlO-8) 

^OO 

0.0 

135 

4.5306 

1.3561 

1.6470 

0.1 

120 

4.5306 

1.3475 

1.6559 

0.3 

124 

4.5301 

1.3574 

1.6479 

0.5 

122 

4.5300 

1.3875 

1.6267 

1.0 

126 

4.5283 

1.4698 

1.5787 

3.0 

118 

4.5234 

1.7950 

1.4231 

5.0 

125 

4.5185 

2.4188 

1.2588 

Table  29:  Parameter  estimation  of  es>o,  To,  and  ks  for  0%,  1%,  3%  and  5%  constant 
variance  noise  (test  2). 


%  RN 

Iter 

re+i  -  j05i 

J6 

HVP'lli. 

FC 

0.0 

135 

7.3419  x  10“9 

112.2835 

0.0281 

244 

0.1 

120 

5.7082  x  10~9 

112.2681 

0.0355 

216 

0.3 

124 

6.8767  x  10~9 

113.8858 

0.0023 

226 

0.5 

122 

2.4453  x  10~9 

117.7016 

0.0287 

217 

1.0 

126 

8.1947  x  10“9 

136.8571 

0.0159 

225 

3.0 

118 

9.0603  x  10~9 

350.8494 

0.0205 

206 

5.0 

125 

4.3675  x  10“9 

784.6225 

0.0296 

231 

coefficient  again  by  adding  constant  variance  noise  to  our  data.  The  other  four 
parameters  are  fixed  at  the  values  given  in  (116),  and  (117).  Table  28  displays  the 
results  of  the  inverse  problem.  Again  we  note  that  the  Nelder-Mead  converges  to  a 
local  minima. Figures  86,  87,  88  and  89  plot  the  details  of  the  Nelder-Mead  algorithm 
and  the  variation  of  the  parameters  with  the  iteration  count.  In  Table  30  we  calculate 
the  confidence  intervals  for  the  three  parameter  estimation  problem.  Again  we  observe 
that  the  intervals  are  larger  as  compared  with  the  results  of  the  Levenberg-Marquardt 
method. 


8  Conclusions  and  future  directions 

In  this  report  we  have  presented  an  electromagnetic  interrogation  technique  for  iden¬ 
tifying  dielectric  parameters  of  a  Debye  medium  by  using  acoustic  pressure  waves 
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Table  30:  Confidence  Intervals  for  the  parameter  esti¬ 
mation  of  600,0,  no  and  for  constant  variance  noise 
with  the  Nelder-Mead  algorithm. 


RN  =  0.0% 

^00, 0 
^0 
^00 

4.5306  ±8.2512  x  10“3 
(1.3561  ±  5.9711  x  10-1)  x  10~8 
1.6470  ±  1.9556  x  10”1 

RN  =  0.1% 

^00, 0 

TO 

^00 

4.5306  ±  8.2283  x  10~3 
(1.3475  ±  5.9323  x  10”1)  x  10“8 
1.6559  ±  1.9523  x  10”1 

RN  =  0.3% 

^00, 0 

T0 

^00 

4.5301  ±  8.2937  x  10“3 
(1.3574  ±6.0170  x  10"1)  x  10~8 
1.6479  ±  1.9682  x  10”1 

RN  =  0.5% 

^00, 0 

T0 

^00 

4.5300  ±  8.4345  x  10“3 
(1.3875  ±  6.2477  x  10”1)  x  10“8 
1.6267  ±2.0063  x  10”1 

RN  =  1.0% 

^00, 0 

T0 

^00 

4.5283  ±  9.0550  x  10“3 
(1.4698  ±  7.1148  x  10"8)  x  10“8 
1.5787  ±2.1736  x  10”1 

RN  =  3.0% 

^00, 0 

T0 

^00 

4.5234  ±  1.4519  x  10~'2 
(1.7950  ±  1.3804)  x  10~8 
1.4231  ±  3.5550  x  10”1 

RN  =  5.0% 

^00, 0 

T0 

^00 

4.5185  ±2.1562  x  10“2 
(2.4188  ±  2.7552)  x  10~8 
1.2588  ±  5.4570  x  10”1 

114 


Function  Differences 


Least  Squares  Function  Value 


0  20  40  60  80  100  120  140 

L°°  Norm  of  Gradient 


Max  Oriented  Lenath 


0  20  40  60  80  100  120  140 


Average  Least  Squares  value 
10  C - T - T - T - T - n 


0  20  40  60  80  100  120  140 


L2  Norm  of  Gradient 

104  | - T - T - T - T - 


0  20  40  60  80  100  120  140 

Iterations 


Figure  86:  Details  of  the  Nelder-Mead  algorithms  for  the  three  parameter  estimation 
°f  £oo,o>  ro  and  Koo  with  varying  levels  (y)  of  constant  variance  noise  added  to  the 
data  as  tabulated  in  Table  28,  and  29. 
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Figure  87:  Variation  in  the  estimate  e^o  over  all  iterations  in  the  Nelder-Mead  algo¬ 
rithm  for  varying  levels  (u)  of  constant  variance  noise  added  to  the  data  as  tabulated 
in  Table  28,  and  29. 


x^=  3.16x  10-8 


Figure  88:  Variation  in  the  estimate  To  over  all  iterations  in  the  Nelder-Mead  algo¬ 
rithm  for  varying  levels  (y)  of  constant  variance  noise  added  to  the  data  as  tabulated 
in  Table  28,  and  29. 
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Kl= 165 


Figure  89:  Variation  in  the  estimate  /too  over  all  iterations  in  the  Nelder-Mead  algo¬ 
rithm  for  varying  levels  (z/)  of  constant  variance  noise  added  to  the  data  as  tabulated 
in  Table  28,  and  29. 


as  virtual  reflectors  for  the  incident  electromagnetic  pulse.  We  considered  the  two 
dimensional  TE  mode  of  Maxwell’s  equations  which  incorporates  the  pressure  de¬ 
pendence  in  the  Debye  medium  via  a  model  for  orientational  polarization.  As  a  first 
approximation  we  assumed  that  the  dielectric  parameters  e^o,  €00,0  and  the  relaxation 
To  are  affine  functions  of  pressure.  We  have  used  the  finite  difference  time  domain 
scheme  to  discretize  our  electromagnetic/acoustic  model  and  to  compute  simulated 
data  to  be  used  in  a  parameter  identification  problem. 

We  formulated  an  inverse  problem  for  the  identification  of  the  material  parameters 
as  well  as  the  coefficients  of  pressure  in  the  Debye  model,  based  on  the  method  of 
least  squares.  We  used  two  different  methods  to  solve  the  inverse  problem;  namely 
Levenberg-Marquardt  and  Nelder-Mead.  The  inverse  problem  indicates  that  the  co¬ 
efficients  that  can  be  identified  depend  on  the  order  of  magnitude  of  ut0.  Thus  in 
the  first  test  problem  we  were  able  to  identify  eSio  and  the  corresponding  pressure 
coefficient  his,  and  in  the  second  test  problem  we  were  able  to  identify  e^o  and  the 
corresponding  pressure  coefficient  It  is  not  possible  to  accurately  identify  all 
the  seven  parameters  that  describe  the  Debye  medium  coupled  with  the  model  for 
pressure  dependence  in  either  of  the  examples  that  we  have  considered  in  this  report. 

We  have  also  computed  confidence  intervals  for  all  the  estimates  obtained  using 
statistical  error  analysis  based  on  the  assumption  of  constant  variance  noise.  We 
showed  that  the  intervals  become  larger  as  the  level  of  added  noise  is  increased. 
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There  are  many  questions  that  still  need  to  be  answered.  We  have  used  a  simple 
linear  model  for  the  pressure  term  in  Maxwell’s  equations.  Also  we  have  assumed  that 
the  effect  of  the  electromagnetic  field  on  the  acoustic  wave  is  negligible.  A  dynamic 
coupled  model  of  the  electromagnetic/acoustic  interaction  is  the  topic  of  a  future 
paper.  We  will  also  consider  the  dielectric  parameters  to  be  nonlinear  functions  of 
pressure  in  a  future  paper  to  determine  if  the  dynamics  of  the  problem  change  with 
the  dependence  of  the  dielectric  parameters  on  the  pressure. 
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